{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple logistic regression\n",
    "\n",
    "This notebook follows John H McDonald's [Handbook of Biological Statistics](http://www.biostathandbook.com/simplelogistic.html) chapter on simple logistic regression.\n",
    "\n",
    "This notebook is provided with a CC-BY-SA license."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import io\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "import scipy.special\n",
    "import seaborn as sns\n",
    "sns.set_style('white')\n",
    "sns.set_context('notebook')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Spider data from *Suzuki et al. (2006)*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Grain size (mm)</th>\n",
       "      <th>Spiders</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.245</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.247</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.285</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.299</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.327</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Grain size (mm) Spiders\n",
       "0            0.245   False\n",
       "1            0.247   False\n",
       "2            0.285    True\n",
       "3            0.299    True\n",
       "4            0.327    True"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = \"\"\"Grain size (mm)\tSpiders\n",
    "0.245\tabsent\n",
    "0.247\tabsent\n",
    "0.285\tpresent\n",
    "0.299\tpresent\n",
    "0.327\tpresent\n",
    "0.347\tpresent\n",
    "0.356\tabsent\n",
    "0.36\tpresent\n",
    "0.363\tabsent\n",
    "0.364\tpresent\n",
    "0.398\tabsent\n",
    "0.4\tpresent\n",
    "0.409\tabsent\n",
    "0.421\tpresent\n",
    "0.432\tabsent\n",
    "0.473\tpresent\n",
    "0.509\tpresent\n",
    "0.529\tpresent\n",
    "0.561\tabsent\n",
    "0.569\tabsent\n",
    "0.594\tpresent\n",
    "0.638\tpresent\n",
    "0.656\tpresent\n",
    "0.816\tpresent\n",
    "0.853\tpresent\n",
    "0.938\tpresent\n",
    "1.036\tpresent\n",
    "1.045\tpresent\n",
    "\"\"\"\n",
    "df = pd.read_table(io.StringIO(data))\n",
    "df.Spiders = df.Spiders == 'present'\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAFmCAYAAACfjbj/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1UVXWi//HP4QAaSsgpQAXrGquH8aFkbI2DUEyCM44r\n762uDwhD6lS3x7GbWj5Niv4gmBl1uuW18tpNTZNrMzJlqxsz0WhdgbQyLDRNKYVQ8YDKg6IC+/eH\n13NFhQMI+3A279dartXe3+13f7bb+rQfOMdmGIYhAABgCT6eDgAAADoOxQ4AgIVQ7AAAWAjFDgCA\nhVDsAABYCMUOAICFeKzY9+3bp9GjR2v9+vWXjRUUFGjSpElKSkrS/PnzPZAOAADv5JFiP336tNLS\n0hQdHX3F8YULF+rll1/WW2+9pZqaGn388ccmJwQAwDt5pNh79OihVatWKTQ09IrjmzZtco05HA6d\nOHHCzHgAAHgtjxS7j4+P/P39mx3v1auXJKm8vFx5eXmKi4szKxoAAF6ty748V1FRoccff1ypqakK\nCgpqdrv6+nqVlpaqvr7exHQAAHRNXbLYa2pq9Mgjj2jGjBnNPoe/4MiRI4qPj9eRI0dMSgcAQNfV\nJYs9MzNT06ZNU0xMjKejAADgVXw9sdOioiJlZmaqrKxMvr6+ysnJ0ahRoxQREaHY2Fi9++67OnTo\nkDZu3CibzaZx48ZpwoQJnogKAIBX8UixDx48WG+++Waz47t27TIxDQAA1tElb8UDAID2odgBALAQ\nih0AAAuh2AEAsBCKHQAAC6HYAQCwEIodAAALodgBALAQih0AAAuh2AEAsBCKHQAAC6HYAQCwEIod\nAAALodgBALAQih0AAAuh2AEAsBCKHQAAC6HYAQCwEIodAAALodgBALAQih0AAAuh2AEAsBCKHQAA\nC6HYAQCwEIodAAALodgBALAQih0AAAuh2AEAsBCKHQAAC6HYAQCwEIodAAALodgBALAQih0AAAuh\n2AEAsBCKHQAAC/FYse/bt0+jR4/W+vXrLxvLy8vThAkTlJiYqBUrVnggHQAA3skjxX769GmlpaUp\nOjr6iuPp6elavny5NmzYoG3btunAgQMmJwQAwDt5pNh79OihVatWKTQ09LKxkpIS9enTR2FhYbLZ\nbIqLi1NBQYEHUgIA4H18PbFTHx8f+fv7X3HM6XTK4XC4lh0Oh0pKSsyK1ikqKo9r8ZI35Kxu0HWB\nPkqd9Ws5HMFux9o7Z2u2u7D+B2eNjh7+Qf0jBios2K/V+29vnqs55s7K0xHzmHFMHZGzM+b19LF7\nO7P/rnqat+RsyaXH8PRDD+il17O7zDF1+ZfnDMPwdISrtnjJGzpQd5Oq/G9RcV2kFi15o1Vj7Z2z\nNdtdWF/8w0kF3fxLnQr4UZv23948bd22vTpqH1f752yWztp/a+b19LF7O7P/rnqat+RsyaXHMPVf\nM7vUMXW5Yg8NDdWxY8dcy0ePHr3iLXtv4qxukM1mkyTZbDY5qxtaNdbeOVuz3YX1vn4927X/9uZp\n67bt1VH7uNo/Z7N01v5bM6+nj93bmf131dO8JWdLLj2Gcz69u9QxdbliDw8PV21trcrKylRfX68t\nW7YoNjbW07GuynWBPq47D4Zh6PpAe6vG2jtna7a7sL7+7Ol27b+9edq6bXt11D6u9s/ZLJ21/9bM\n6+lj93Zm/131NG/J2ZJLj8GvobpLHZM9NTU11eydFhUVaebMmdqxY4e++uor/e1vf9PJkyfldDp1\n00036ZZbbtHChQu1adMmjR07Vj/72c+anauqqkpr167VlClTdO2115p3EG0wIupH2rX972o47VRE\nr5NaOGuarrnmGrdj7Z2zNdtdWO/r66OKg1/IcU29bgisavX+25vnao65s/J0xDxmHFNH5OyMeT19\n7N7O7L+rnuYtOVty6TFkzH9E3+7K7zLHZDO8/CF2aWmp4uPjlZubq4iICE/HAQDAo7rcrXgAANB+\nFDsAABZCsQMAYCEUOwAAFkKxAwBgIRQ7AAAWQrEDAGAhFDsAABZCsQMAYCEUOwAAFkKxAwBgIRQ7\nAAAWQrEDAGAhFDsAABZCsQMAYCEUOwAAFkKxAwBgIRQ7AAAWQrEDAGAhFDsAABZCsQMAYCEUOwAA\nFkKxAwBgIRQ7AAAWQrEDAGAhFDsAABZCsQMAYCEUOwAAFkKxAwBgIRQ7AAAWQrEDAGAhFDsAABZC\nsQMAYCEUOwAAFkKxAwBgIb6e2GlGRoYKCwtls9k0b948DR061DW2fv16bd68WXa7XUOGDNHcuXM9\nEREAAK9kerHv2LFDBw8eVFZWlg4cOKD58+crKytLklRTU6PXX39dubm5stlseuihh7Rr1y7dfvvt\nZscEAMArmX4rPj8/XwkJCZKkyMhIVVVVqba2VpLk7+8vf39/1dTUqL6+XnV1dQoKCjI7IgAAXsv0\nYnc6nXI4HK7l4OBgOZ1OSeeL/cknn1RCQoLi4+N1++2368YbbzQ7IgAAXsvjL88ZhuH655qaGr32\n2mv661//qtzcXBUWFmrv3r0eTAcAgHcxvdhDQ0NdV+iSVF5erpCQEElScXGxBgwYoKCgIPn6+urO\nO+9UUVGR2REBAPBaphd7TEyMcnJyJElFRUUKCwtTQECAJCk8PFzFxcU6e/asJOnrr7/mVjwAAG1g\n+lvxUVFRGjx4sBITE2W327VgwQJlZ2crMDBQCQkJeuihh5SSkiJfX19FRUVp+PDhZkcEAMBr2YyL\nH3J7odLSUsXHxys3N1cRERGejgMAgEd5/OU5AADQcSh2AAAshGIHAMBCKHYAACyEYgcAwEIodgAA\nLIRiBwDAQih2AAAshGIHAMBCKHYAACyEYgcAwEIodgAALIRiBwDAQih2AAAshGIHAMBCKHYAACyE\nYgcAwEIodgAALIRiBwDAQih2AAAsxG2xnzhxwowcAACgAzRb7Nu2bVN0dLRGjx6tSZMmqbKy0sxc\nAACgHZot9uXLl2vjxo3asWOH/vEf/1FPPvmkzp49a2Y2AADQRs0Wu5+fnwYMGCBJSk5OVkJCgsaO\nHavS0lKlpaWZFhAAALRes8VuGIaKi4tdyw899JA2btyokJAQjR8/3pRwAACgbZot9meffVYPP/yw\n8vLyXOscDod69Oih2267zZRwAACgbZot9ttvv125ubkKCwu7bOzLL7/s1FAAAKB9mi32qqoqlZaW\nat68eSopKXH9Ki4u1uzZs83MCAAAWsm3uYGdO3dqzZo12rNnj6ZMmeJa7+Pjo9jYWFPCAQCAtmm2\n2OPi4hQXF6cNGzZo8uTJZmYCAADt1GyxX5CQkKA1a9bo5MmTMgzDtf7pp5/u1GAAAKDt3H6k7KOP\nPqpvvvlGPj4+stvtrl8AAKDrcXvFHhAQoIyMDDOyAACAq+T2iv2OO+7QgQMHzMgCAACuktsr9k8+\n+USrV69WcHCwfH19ZRiGbDabtmzZYkI8AADQFm6L/ZVXXunwnWZkZKiwsFA2m03z5s3T0KFDXWNH\njhzRjBkzVF9fr0GDBik1NbXD9w8AgFW5vRUfEhKiLVu2aMOGDQoPD5fT6dT111/f7h3u2LFDBw8e\nVFZWltLS0pSent5kPDMz0/W59Ha7XUeOHGn3vgAA6G7cFntqaqoOHTqkTz/9VJJUVFSkOXPmtHuH\n+fn5SkhIkCRFRkaqqqpKtbW1ks5/8cznn3+uUaNGSZKef/559e3bt937AgCgu3Fb7MXFxZo7d656\n9uwpSUpKSlJ5eXm7d+h0OuVwOFzLwcHBcjqdkqTKykoFBAQoPT1dSUlJWrZsWbv3AwBAd+S22H19\nzz+Gt9lskqRTp06prq6uwwJc/KE3hmGovLxcU6dO1bp167R7925t3bq1w/YFAIDVuS32MWPGaMqU\nKSotLVVaWpruu+8+jRs3rt07DA0NdV2hS1J5eblCQkIknb96Dw8PV0REhHx8fBQdHa39+/e3e18A\nAHQ3bov9V7/6lWbOnKmkpCTdcMMNWrZsmaZOndruHcbExCgnJ0fS+ef1YWFhCggIkCTZ7XZFRETo\n0KFDrvGBAwe2e18AAHQ3bn/c7eTJk7rmmmv08MMP6+OPP9bWrVsVFhbmuspuq6ioKA0ePFiJiYmy\n2+1asGCBsrOzFRgYqISEBM2bN09z5syRYRi65ZZbXC/SAQAA92zGxQ+5r+Bf/uVfNGXKFIWHh+up\np57S5MmTtXXrVq1cudKsjC0qLS1VfHy8cnNzFRER4ek4AAB4lNtb8adPn1ZMTIw++OADJScnKzk5\nWefOnTMjGwAAaKNWFXtlZaVycnL0s5/9TIZh6OTJk2ZkAwAAbeS22MeNG6ef//zn+ulPf6p+/frp\n3//93zVixAgzsgEAgDZy+4z9UlVVVbr22ms7K0+b8YwdAID/4/aK/ZtvvtEDDzygMWPGSJLefPNN\nFRYWdnowAADQdm6LffHixXrhhRdcP942duxYZWRkdHowAADQdq36SNnbbrvNtTxw4EDXx8wCAICu\npVXFXlJS4vqs+K1bt6qNj+UBAIBJ3F56z549W0888YS+++47DR8+XOHh4fr9739vRjYAANBGbos9\nODhYmzdvVmVlpfz9/dW7d28zcgEAgHZweyt+1qxZkiSHw0GpAwDQxbm9Yv+Hf/gHPffcc4qKipKf\nn59r/fjx4zs1GAAAaDu3xX7u3DnZ7Xbt2rWryXqKHQCArsdtsV/4mfWKigrZbDY5HI5ODwUAANrH\nbbG///77Sk9Pl81mk2EYru9QT0hIMCMfAABoA7fF/uqrr2rDhg264YYbJEnfffednn76aYodAIAu\nyO1b8SEhIa5Sl85/8hxftgIAQNfk9or95ptvVlpamu666y41NjaqoKBA/fr1U35+viQpOjq600MC\nAIDWcVvsRUVFkqS9e/c2Wb9v3z7ZbDaKHQCALsRtsb/55ptm5AAAAB3A7TN2AADgPSh2AAAspM3F\n3tjY2Bk5AABAB3Bb7Js2bdL69etVX1+vyZMnKz4+Xm+99ZYZ2QAAQBu5Lfb/+q//0oQJE/Thhx/q\n5ptvVm5urv77v//bjGwAAKCN3BZ7jx495O/vr61bt+qXv/ylfHx4LA8AQFfVqpZetGiRvvjiC/3k\nJz/Rzp07dfbs2c7OBQAA2sFtsS9ZskQ33nijXnnlFdntdv3www9atGiRGdkAAEAbuf2Amv/4j//Q\n/PnzXcv33ntvpwYCAADt5/aK3W63Kz8/X2fOnFFjY6PrFwAA6HrcXrG//fbbWrNmjQzDcH0nu81m\n0549e8zIBwAA2sBtsX/++edm5AAAAB3A7a34kydP6ne/+52effZZSdJHH32kysrKTg8GAADazm2x\n//a3v1W/fv1UUlIiSTp79qxmz57d6cEAAEDbuS32yspKPfjgg/Lz85MkjRkzRnV1dZ0eDAAAtF2r\nPqDm3LlzstlskiSn06lTp05d1U4zMjKUmJioyZMn66uvvrriNkuXLlVKSspV7QcAgO7G7ctzycnJ\nGj9+vI4dO6bHHntMX331VZOfa2+rHTt26ODBg8rKytKBAwc0f/58ZWVlNdnmwIED+uyzz1x3CQAA\nQOu4LfaxY8fqxz/+sXbu3Cl/f38tXrxYoaGh7d5hfn6+EhISJEmRkZGqqqpSbW2tevXq5domMzNT\nM2bM0Msvv9zu/QAA0B01W+x/+ctfLlt35swZ5eXlSZLuu+++du3Q6XRqyJAhruXg4GA5nU5XsWdn\nZ2vEiBHq379/u+YHAKA7a7bYt23bJkk6fvy4vvnmG91xxx1qaGjQrl27FBUV1e5iv5RhGK5/Pnny\npDZt2qTVq1fr8OHDTcYAAIB7zRb7H/7wB0nS9OnT9eGHH6pnz56SpJqaGv32t79t9w5DQ0PldDpd\ny+Xl5QoJCZEkFRQU6Pjx40pOTtaZM2dUUlKizMxMzZkzp937AwCgO3H7VnxZWZmr1CWpd+/eKisr\na/cOY2JilJOTI0kqKipSWFiYAgICJEm/+MUv9N577ykrK0vLly/XoEGDKHUAANrA7ctzN998sxIT\nExUVFSUfHx8VFhbqhhtuaPcOo6KiNHjwYCUmJsput2vBggXKzs5WYGCg66U6AADQPjbDzYNswzCU\nl5enffv2yTAMRUZG6q677pKPT6t+BL7TlZaWKj4+Xrm5uYqIiPB0HAAAPKrZdt69e7ek88+9fXx8\ndNttt+lHP/qR/P399emnn5oWEAAAtF6zt+LfeecdDRo0SCtWrLhszGazKTo6ulODAQCAtnN7K76r\n41Y8AAD/x+2D8u3bt+uBBx7QHXfcoWHDhmnSpEn68ssvzcgGAADayO1b8S+88IJmz56t4cOHyzAM\nffbZZ0pNTb3iJ9MBAADPcnvF3qdPH0VHR8vf3189evRQTEyMwsLCzMgGAADayO0V+x133KHVq1cr\nNjZWjY2NKigoUGRkpEpKSiRJAwYM6PSQAACgddwW++bNmyVJa9eubbL+gw8+kM1mU25ubuckAwAA\nbea22D/66CMzcgAAgA7Q7DP2mpoarV692rWclZWlf/qnf9L06dObfIkLAADoOpot9gULFqiiokKS\n9N1332nZsmWaPXu2Ro4cqfT0dNMCAgCA1mu22EtKSjRz5kxJUk5OjsaMGaORI0cqMTGRK3YAALqo\nZov9wlepSuc/pOanP/2pa9lms3VuKgAA0C7NFntDQ4MqKip06NAh7dy5UzExMZKk2tpanT592rSA\nAACg9Zp9K/6RRx7R2LFjVVdXp6eeekpBQUGqq6tTUlKSJk6caGZGAADQSs0We1xcnP7nf/5HZ86c\nUe/evSVJPXv21LPPPqvY2FjTAgIAgNZr8efY/fz85Ofn12QdpQ4AQNfl9rPiAQCA96DYAQCwEIod\nAAALodgBALAQih0AAAuh2AEAsBCKHQAAC6HYAQCwEIodAAALodgBALAQih0AAAuh2AEAsBCKHQAA\nC6HYAQCwEIodAAALodgBALAQih0AAAvx9cROMzIyVFhYKJvNpnnz5mno0KGusYKCAv3xj3+U3W7X\nwIEDlZ6e7omIAAB4JdOv2Hfs2KGDBw8qKytLaWlplxX3woUL9fLLL+utt95STU2NPv74Y7MjAgDg\ntUwv9vz8fCUkJEiSIiMjVVVVpdraWtf4pk2bFBoaKklyOBw6ceKE2REBAPBaphe70+mUw+FwLQcH\nB8vpdLqWe/XqJUkqLy9XXl6e4uLizI4IAIDX8vjLc4ZhXLauoqJCjz/+uFJTUxUUFOSBVAAAeCfT\niz00NLTJFXp5eblCQkJcyzU1NXrkkUc0Y8YMRUdHmx0PAACvZnqxx8TEKCcnR5JUVFSksLAwBQQE\nuMYzMzM1bdo0xcTEmB0NAACvZ/qPu0VFRWnw4MFKTEyU3W7XggULlJ2drcDAQMXGxurdd9/VoUOH\ntHHjRtlsNo0bN04TJkwwOyYAAF7JZlzpIbcXKS0tVXx8vHJzcxUREeHpOAAAeJTHX54DAAAdh2IH\nAMBCKHYAACyEYgcAwEIodgAALIRiBwDAQih2AAAshGIHAMBCKHYAACyEYgcAwEIodgAALIRiBwDA\nQih2AAAshGIHAMBCKHYAACyEYgcAwEIodgAALIRiBwDAQih2AAAshGIHAMBCKHYAACyEYgcAwEIo\ndgAALIRiBwDAQih2AAAshGIHAMBCKHYAACyEYgcAwEIodgAALIRiBwDAQih2AAAshGIHAMBCKHYA\nACyEYgcAwEIodgAALMTXEzvNyMhQYWGhbDab5s2bp6FDh7rG8vLy9Mc//lF2u1133323nnjiCU9E\nBADAK5l+xb5jxw4dPHhQWVlZSktLU3p6epPx9PR0LV++XBs2bNC2bdt04MABsyMCAOC1TC/2/Px8\nJSQkSJIiIyNVVVWl2tpaSVJJSYn69OmjsLAw2Ww2xcXFqaCgwOyIAAB4LdNvxTudTg0ZMsS1HBwc\nLKfTqV69esnpdMrhcLjGHA6HSkpKzI6oisrjmvf/VujLPQflF3Cthtx0nV6Y95gcjmDTs3S0isrj\nWrzkDTmrG3RdoI+efugBvfR6tms5ddav23ycl855YY7m1rdlju6iM84LgO7J4y/PGYbRrrHOtHjJ\nG9q+t1Ihg8YqeOBdKm38kRYtecMjWTra4iVv6EDdTaryv0XFdZGa+q+ZTZbbc5yXznlhjubWt2WO\n7qIzzguA7sn0Yg8NDZXT6XQtl5eXKyQkxDV27Ngx19jRo0cVGhpqdkQ5qxvk599TNptNkmSz2eSs\nbjA9R2dwVjc0Oa5zPr2v+jgvnfPCHM2tb8sc3UVnnBcA3ZPpxR4TE6OcnBxJUlFRkcLCwhQQECBJ\nCg8PV21trcrKylRfX68tW7YoNjbW7Ii6LtBH586edt0xMAxD1wfaTc/RGa4L9GlyXH4N1Vd9nJfO\neWGO5ta3ZY7uojPOC4DuyWZ44H73smXLtH37dtntdi1YsEC7d+9WYGCgEhIS9Nlnn2nJkiWSpDFj\nxmjq1KktzlVaWqr4+Hjl5uYqIiKiQ/JVVh7X3P99xu4fcK2G3HS90uc9aolnnJWVx7Xof5/lXh9o\n1/SH79dLq7JdywtnTWvzcV4654U5mlvfljm6i844LwC6J48Ue0fqjGIHAMBbefzlOQAA0HEodgAA\nLIRiBwDAQih2AAAshGIHAMBCKHYAACyEYgcAwEIodgAALIRiBwDAQih2AAAshGIHAMBCKHYAACyE\nYgcAwEIodgAALIRiBwDAQih2AAAshGIHAMBCKHYAACyEYgcAwEIodgAALIRiBwDAQih2AAAshGIH\nAMBCKHYAACyEYgcAwEIodgAALIRiBwDAQih2AAAshGIHAMBCKHYAACyEYgcAwEIodgAALIRiBwDA\nQih2AAAshGIHAMBCfM3eYX19vebMmaOysjLZ7XZlZGQoIiKiyTbvv/++3njjDdntdo0YMULPPPOM\n2TEBAPBKpl+xv/feewoKCtJbb72lxx57TEuXLm0yXldXp6VLl2rt2rXKyspSfn6+Dhw4YHZMAAC8\nkunFnp+fr4SEBEnSyJEj9cUXXzQZ79mzpzZv3qxrrrlGktSnTx+dOHHC7JgAAHgl04vd6XTK4XBI\nkmw2m3x8fFRfX99km4CAAEnS3r17VVZWpmHDhpkdEwAAr9Spz9jffvtt/elPf5LNZpMkGYahXbt2\nNdmmsbHxir/3+++/16xZs7R06VLZ7fZm99HQ0CBJOnLkSAelBgCg6+vbt698fS+v8U4t9gkTJmjC\nhAlN1s2dO1dOp1O33nqr60r90mBHjhzRb37zG/3hD3/Qrbfe2uI+jh07JklKTk7uwOQAAHRtubm5\nl718LnngrfiYmBh98MEHiomJ0UcffaQRI0Zcts38+fO1cOFC3XbbbW7nGzJkiNavX6+QkJAWr+wB\nALCSvn37XnG9zTAMw8wgjY2Nmj9/vg4ePKgePXooMzNTYWFhWrlypUaMGKGgoCDdf//9Gjp0qAzD\nkM1m07Rp03TPPfeYGRMAAK9kerEDAIDOwyfPAQBgIRQ7AAAWQrEDAGAhXlXsGRkZSkxM1OTJk/XV\nV181GSsoKNCkSZOUlJSk+fPneyhh99LS+bhg6dKlSklJMTlZ99PSuThy5IiSkpI0ceJEpaameiZg\nN9LSuVi/fr0SExOVnJysjIwMDyXsXvbt26fRo0dr/fr1l43l5eVpwoQJSkxM1IoVKzyQrpMYXmL7\n9u3Go48+ahiGYezfv9+YNGlSk/Gf//znxtGjRw3DMIzp06cbW7duNT1jd+LufFxYn5iYaKSkpJgd\nr1txdy6efvpp48MPPzQMwzAWL15sHD582PSM3UVL56K6utq45557jMbGRsMwDOPXv/61UVhY6JGc\n3cWpU6eMlJQU4/nnnzfWrVt32fjYsWONI0eOGI2NjUZSUpKxf/9+D6TseF5zxX7xZ8xHRkaqqqpK\ntbW1rvFNmzYpNDRUkuRwOPh8+U7m7nxIUmZmpmbMmOGJeN1KS+fCMAx9/vnnGjVqlCTp+eefb/Zn\nX3H1WjoX/v7+8vf3V01Njerr61VXV6egoCBPxrW8Hj16aNWqVa5uuFhJSYn69OmjsLAw2Ww2xcXF\nqaCgwAMpO57XFPvFnzEvScHBwXI6na7lXr16SZLKy8uVl5enuLg40zN2J+7OR3Z2tkaMGKH+/ft7\nIl630tK5qKysVEBAgNLT05WUlKRly5Z5Kma30NK58Pf315NPPqmEhATFx8fr9ttv14033uipqN2C\nj4+P/P39rzh26blyOBwqLy83K1qn8ppiv5RxhR+/r6io0OOPP67U1FT+T9hkF5+PkydPatOmTZo2\nbZoMw7jiuULnufjP2zAMlZeXa+rUqVq3bp12796trVu3ejBd93LxuaipqdFrr72mv/71r8rNzVVh\nYaH27t3rwXS4mJX+O+U1xR4aGtrkirC8vFwhISGu5ZqaGj3yyCOaMWOGoqOjPRGxW2npfBQUFOj4\n8eNKTk7Wb37zG+3Zs0eZmZmeimp5LZ2L4OBghYeHKyIiQj4+PoqOjtb+/fs9FdXyWjoXxcXFGjBg\ngIKCguTr66s777xTRUVFnora7YWGhrq+a0SSjh49esVb9t7Ia4o9JiZGOTk5kqSioiKFhYW5vt5V\nOv88d9q0aYqJifFUxG6lpfPxi1/8Qu+9956ysrK0fPlyDRo0SHPmzPFkXEtr6VzY7XZFRETo0KFD\nrvGBAwd6LKvVtXQuwsPDVVxcrLNnz0qSvv76a27Fe1B4eLhqa2tVVlam+vp6bdmyRbGxsZ6O1SG8\n6iNlly1bpu3bt8tut2vBggXavXu3AgMDFRsbq5/85CcaNmyY6/Plx40bd9k3y6FjNXc+Lrw8JEk/\n/PCD5s6dq7Vr13owqfW1dC4OHTqkOXPmyDAM3XLLLVq0aJGn41paS+di48aN+vOf/yxfX19FRUVp\n1qxZno5raUVFRcrMzFRZWZl8fX0VFhamUaNGKSIiQgkJCfrss8+0ZMkSSdKYMWM0depUzwbuIF5V\n7AAAoGVecyseAAC4R7EDAGAhFDsAABZCsQMAYCEUOwAAFkKxAwBgIRQ74AWcTqfmzJmj+++/X8nJ\nybr//vsWDeyHAAAFRElEQVTb9dkAGRkZ2r1791VlWblyZYd8LG15ebkSExN18uTJq57rgunTp2vb\ntm0dNh/gjfg5dsALTJw4UePHj9fEiRMlnf9ehKlTp2r69OkaPXq0h9O1z6OPPqopU6Zo5MiRHTbn\niRMnNHHiRL3zzju65pprOmxewJv4ejoAgJbl5+fL19fXVeqSdN111yk7O1u+vuf/FZ47d678/Pz0\n/fffa8mSJdq1a5dWrVqlHj16qKGhQb///e/Vv39/paSk6IknnpDdbtfKlSvVt29f7d+/X35+fq7t\nLzh16pRmzpyp6upq1dfX65577tGjjz6quXPnavjw4XI4HHrjjTdks9l06tQpffvttyosLFRVVZUW\nLlyo48ePq7q6WtOmTdO9997b5Jj27Nmjw4cPu0p91KhRmjx5sj755BMdO3ZMs2fPVlZWloqLi/XE\nE0/ovvvu09y5c9WnTx8VFxdr//79mjlzpj766CPt3btXw4cPV2pqqvr06aN77rlHb7/9th588EET\nzg7Q9XArHujivv32Ww0ZMuSy9RdK/YK6ujqtXbtWoaGhqq6u1osvvqg1a9bo7rvv1rp16y77/YWF\nhZo5c6aysrJks9n0ySefNBnPy8tTQ0OD1q1bpw0bNiggIKDJN2CNGjVKb775ptauXauIiAg9//zz\nkqQXX3xRd999t1avXq1169bppZde0vHjx5vM/cknn+iuu+5qsu66667T2rVrNWzYMK1du1avvvqq\n0tLStGbNGtc2lZWVeu211/TUU09p8eLFSk1N1Z/+9CdlZ2erpqZGkjRy5MjLjgXoTrhiB7o4u92u\nhoYG1/LGjRu1efNmnT17Vv369dOLL74oSYqKinJt43A49Nxzz8kwDDmdTg0bNuyyeSMjIxUcHCzp\n/BdiXPqs+8c//rFeeuklPfPMM7r77rs1YcIE2Wy2y+Z5/fXXFRgYqPHjx0uSPv30U3399dfatGmT\npPPfQ15aWuralyQdPnxYkZGRTea5kD8sLEx9+/aVJPXt21fV1dVNMl1YHxkZqd69e0s6/y121dXV\n6t27t/r3768ffvih+T9QwOIodqCLu/XWW/XnP//ZtTxx4kRNnDhR27dv17/927+51vv5+UmS6uvr\n9cwzz+idd97RgAEDtH79en399deXzWu325ssX/q6jcPh0LvvvqudO3fqww8/1AMPPKC//OUvTbYp\nKCjQ3/72tyZ3BPz9/bVw4UINHjy4Tcd58R2Ii7NdnOvi9Zduw+tCwHncige6uDvvvFPBwcFauXKl\na925c+e0bds29ezZ87Lta2trZbfb1b9/f505c0a5ubmurwpti23btunvf/+7oqKi9Oyzz6pXr16q\nqKhwjR89elRpaWl68cUXm5Ty8OHD9f7770s6/3hg0aJFamxsbDJ3v379dPjw4TZnao2ysjKFh4d3\nytyAN6DYAS/wyiuvqKKiQvfdd59SUlKUmJiouro611dOXiwoKEj33nuv/vmf/1kzZszQww8/rE8/\n/VQ5OTlXvJXenIEDB+o///M/9atf/UoPPvigYmNj1a9fP9f4ihUrVFtbq+eee04pKSl68MEH9f33\n3+upp57SwYMHlZSUpJSUFA0aNEg+Pk3/U3PXXXc1eQ7ellxXcvHvz8vLu+z5PdCd8ONuADziscce\nU0pKimJiYjpszuPHjysxMVHZ2dkKCAjosHkBb0KxA/CIY8eOafr06Xr11VcVFBTUIXNOnz5dkydP\nVnR0dIfMB3gjih0AAAvhGTsAABZCsQMAYCEUOwAAFkKxAwBgIRQ7AAAWQrEDAGAh/x938U1uVzsH\nWgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f884cf60f60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.plot.scatter('Grain size (mm)', 'Spiders')\n",
    "plt.ylabel('Spiders present?')\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I will analyse this with the *scikit-learn* package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sklearn.linear_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "scikit-learn has a logisitic regression classifier which uses regularization. To eliminate regularization, we set the regularization parameter `C` to $10^{12}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.64748882] [[ 5.12118846]]\n"
     ]
    }
   ],
   "source": [
    "# C=1e12 is effectively no regularization - see https://github.com/scikit-learn/scikit-learn/issues/6738\n",
    "clf = sklearn.linear_model.LogisticRegression(C=1e12, random_state=0)\n",
    "clf.fit(df['Grain size (mm)'].reshape(-1, 1), df['Spiders'])\n",
    "print(clf.intercept_, clf.coef_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is in agreement with the equation John reports:\n",
    "$$\n",
    "probability of spider presence = \\frac{e^{-1.6476+5.1215(grain \\; size)}}{(1+e^{-1.6476+5.1215(grain \\; size)}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_log_reg(x, y, data, clf, xmin=None, xmax=None, alpha=1, ax=None):\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots()\n",
    "    else:\n",
    "        fig = ax.figure\n",
    "    ax.scatter(data[x], data[y], color='black', zorder=20, alpha=alpha)\n",
    "    if xmin is None:\n",
    "        xmin = x.min()\n",
    "    if xmax is None:\n",
    "        xmax = x.max()\n",
    "    X_test = np.linspace(xmin, xmax, 300)\n",
    "\n",
    "    loss = scipy.special.expit(X_test * clf.coef_ + clf.intercept_).ravel()\n",
    "    ax.plot(X_test, loss, linewidth=3)\n",
    "\n",
    "    ax.set_xlabel(x)\n",
    "    ax.set_ylabel(y)\n",
    "    fig.tight_layout()\n",
    "    sns.despine()\n",
    "    return fig, ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGBCAYAAACTuDAhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9//H3zGQjIYYEkgAJArIpixBAMARBINSlWpcK\nRCwo7lu1YlUWWeQrhda11p/WtYqggApuVWmJgpYEWQUJSNhJCJCVJXtm5v7+iAxEICSQzJ2bvJ6P\nRx6Te+7Nnc/Jcuede++cYzMMwxAAAIBF2c0uAAAA4FwQZgAAgKURZgAAgKURZgAAgKURZgAAgKUR\nZgAAgKWZFmbS09M1fPhwzZs376R1K1eu1KhRozR69GhNnjzZhOoAAIBVmBJmSkpK9PTTTys+Pv6U\n66dNm6Z//OMfev/991VYWKjvvvvOyxUCAACrMCXMBAYG6s0331RUVNQp1y9atMizLiIiQocOHfJm\neQAAwEJMCTN2u10BAQGnXR8SEiJJys7OVkpKigYPHuyt0gAAgMX47A3AeXl5uu+++zR9+nSFhYWd\ndjun06nMzEw5nU4vVgcAAHyFT4aZwsJC3XXXXRo/fvxp76s55sCBAxo2bJgOHDjgpeoAAIAv8ckw\nM3v2bI0bN04JCQlmlwIAAHycnxlPmpaWptmzZysrK0t+fn5asmSJhg4dqtjYWA0cOFCfffaZ9u7d\nq4ULF8pms+naa6/ViBEjzCgVAAD4OFPCTLdu3fTee++ddv3GjRu9WA0AALAyn7zMBAAAUFOEGQAA\nYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmE\nGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAA\nYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmE\nGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmEGQAAYGmmhZn0\n9HQNHz5c8+bNO2ldSkqKRowYoaSkJL3yyismVAcAAKzClDBTUlKip59+WvHx8adcP3PmTL388sv6\n4IMPtGLFCu3YscPLFQIAAKswJcwEBgbqzTffVFRU1EnrMjIy1KxZM0VHR8tms2nw4MFauXKlCVUC\nAAArMCXM2O12BQQEnHJdbm6uIiIiPMsRERHKzs72VmloIPLy8jRq1Cj169dPo0aN0vbt26ss5+fn\nV7t9fn7+Kfdx3XXXKSIiQhEREbr++utP2k91+6uuvbb98dbznm0d1W1fV7V407l+/71dR10+V1xc\nnNq0aaPevXtb5ueF49xuQ06XW2UVLpWUOVVYUqEjReUqOFqqvMMlyi4o1oG8ImXlFCrj4FHt2X9E\nO/cd1raMAm3dk6+te/JV4XSb3Y0z8jO7gDMxDMPsEmBB999/vxYuXChJWr16tVJSUpSZmelZlqQF\nCxacdvtjTrcPSfr000913333VdlPdftbsGDBadtr25/TfV1dP+/Z1lHd9pLqpBZvOtfvv7frqOvn\nkqTMzEytX7++Xp/TlxmGIZfbUHmFS+UVbpU7Xapwun9Zdqnc6VZFRWVocDrdcrrdlY8ut5wu45fH\nEz6vss3x9a5fPq9wueU6tv0J27oNQy5XZS3HPvc8ut1V29yG6uIlNDK8if75xDAF+DvOfWf1xOfC\nTFRUlHJycjzLBw8ePOXlKKA6u3btqrL86/8mf73+TMun2sfptqtufzV5ntrsr76f92zrqM32Z1uL\nN53r99/bddTHc3njOc+WYVQGgJIyl0rLnSorrzwLUVbuUkm5U2VllY+l5U6V/rJNablLpSdsU/FL\nQDkxmFRddsndSP+3zikoUc6hEsVENjW7lNPyuTATExOjoqIiZWVlKSoqSsuWLdNzzz1ndlmwmPbt\n21c5ExAREaHi4uIq66vb/tj66vZxqv2caX+na69tf7z1vGdbx5m2r4tavOlcv//erqM+nqs+n7PC\n6dKRonIVllSouMSpotIKFZdWqKjUqeKSChWVVqiopELFpcfWOX9ZrtympMwpd2NNGjVgt9vk+OXD\nftKj/fjnNpscjqqPfg67BlzcyqeDjGRSmElLS9Ps2bOVlZUlPz8/LVmyREOHDlVsbKwSExM1bdo0\njR8/XpJ0zTXXqG3btmaUCQt79dVXJVX+F9m+fXv95S9/0aRJkzzLx9afbvsT15+4j0cffVTfffed\nJGnw4MEn7edM+6vueWrTH28979nWUZPtz7UWbzrX77+366jL50pPT1dubq4iIyPVqVOnMz5neYVL\nh46W6VBhmY4UletocbmOFpXryC+PR4srji//0lZa7qq3ftQlu92mQH+7/P0cCvCzy9/foUB/h/z9\n7Ar45THQ3yE/h13+fnY5HJVhwN9hl8Nhl5/DJj8/u/wcxz5sJ3/uZ5fDXvn1fg6bHL98vZ+jcn/H\ng4j9NOHk+Hq73Sa7TbLZbGZ/6+qdzbD4TSmZmZkaNmyYkpOTFRsba3Y5ANDguN2GDhWWKe9wiQ4d\nLVPB0bJfHkurLB86WqqiUqcpNTrsNgUF+qlJgEOBAX4KCnQoKMBPQQEOBQVWPjYJ8FOgZ7nqugD/\nyoAS4H/8c38/hwL87Z5lh4NxZn2Vz11mAgB4j9tt6HBhmXIOlSjvcOW9EXmHSpV7qMTTlne4VK56\nvIzjsNsUGhKgpk38FRLkr+AgPwWf8HlIk8rHpk38FRz063Z/NQn0k78fQaMxI8wAQANXVFKhA3lF\nOpBfrAO5vzzmFelgXrFyDhXL6aq7oOKw2xTWNFDNQgN1XkiAzgsOUGhIgEKDAxQa4l9l+byQyo8m\ngX6N4lII6g9hBgAagNJyp/ZlFyozu1AZ2UeVlVNUGWDyinW0uPyc9x8aHKDmYUGKOC9IzUIDFR5a\nGViahQYpvGmgmp0XqGZNAxUaHCC7nWAC7yLMAICFHC4sqwwsB496gktmdqFyCorPekyRpk381aJZ\nkxM+gtQirPLzyGZNFBEWpKAAXi7gu/jtBAAf5HK5tS+nULuyjmhX1mHt2n9Eu/YdVsHRslrvy9/P\nrpbNgxUdEaKWzYPVsnmIWkYEq2WLEEWFB6tJIC8FsDZ+gwHAZBVOt3bvP6z0vYe0I/OQdu0/or37\nj6i8FsPI2+02tWoerNioUMVGNVVsVFO1atFULZsHKzw0iEs/aNAIMwDgRYZhaH9ekdL3HlL63gKl\n7y3Qzn2Hazz/TYC/Q22im6rNsdASXfnYukWI/P18d7h5oD4RZgCgHlU4XUrfe0hpO/OUtitP2/YW\n6GhxRY2+tnlYkNq3DlP71ud5Hlu1aCoHZ1mAKggzAFCHSsqc+nl3vtJ25mnTzjyl7y2o0VmX6Ihg\ndT4/XJ3aNNMFrcPUrvV5Cmsa6IWKAesjzADAOXC53NqWcUjr03P0Y3q2tu4pOOMAc02b+Kvz+eG/\nfDRTpzbhahZKcAHOFmEGAGopK6dQ67dma316jn7akaviMwzh36pFiLpf0Fxd2zdX1/YRatUihEHi\ngDpEmAGAM3C53Nq8O1+r0g5o9eYD2pdTVO327Vqdp24XNPd8RJwX5KVKgcaJMAMAp1BYUqF1Px/U\nqrSDWvvzQRWWnP6m3RZhQerVOUq9OkeqZ6dILhkBXkaYAYBfFBaXa+WmA1qxMUs/pmefds6iwACH\nenaMVK/OlR+xUU25bASYiDADoFErKqnQD2n79f2P1QeYFmFBuqRrS/Xr1lIXd2yhAH/GdAF8BWEG\nQKNT4XRpVdpBfbs2Q2t/zpbTdeq3TneMDVO/bq3Ur2u0LogJ4+wL4KMIMwAaBcMwtC3jkJau3qvv\n1+877T0wHds002U9WyuhZ4yiI4K9XCWAs0GYAdCg5R0u0TdrMvTNmgxlZheecpuOsWEa2DNGCT1b\nq2XzEC9XCOBcEWYANDhut6Eft+XoyxW7tHrzAZ1qDLuoiGAN7dNGQ/rGqnWLpt4vEkCdIcwAaDCO\nFpcrefVefZmyW/tzTx4LJijAoYSerTXskvPVrX1zZpIGGgjCDADL25ZRoC9X7NZ36zNVfop5kC7u\n2ELDLjlfA3q0UlAghz2goeGvGoAlud2GVm0+oMXLtmvzrvyT1ocE+WlYv/N1VXw7xUaFmlAhAG8h\nzACwlNJyp75dk6FPlu9Q1ikuJV0QE6bfJrTXoF4xnIUBGgn+0gFYQmFxuT7/3y598b+dOlJUXmWd\nw27TZXExuiahvTqfH854MEAjQ5gB4NMOF5bp0+926N8rdp00O3VIkJ+ujG+nay+7QM3DmphUIQCz\nEWYA+KSCo6VavGyHvkrZpdJyV5V1keFNdN2gDhre73wFB/mbVCEAX0GYAeBTjhSV6+NvtumLFbtU\nXlE1xMRENtXIxM4aHBcjh8NuUoUAfA1hBoBPKClz6tPvdmjxsu0nXU5q2zJUo4Z30YCLW8vB2DAA\nfoUwA8BUFU6XvkrZrYXJ6TpcWPXG3gtiwpQ0vLP6d2vFAHcAToswA8AUhmHofz9m6Z1/pym7oKTK\nupjIphpz1UUacHEr3pkE4IwIMwC8Ln1vgd78dJO27K462F1keBON/k0XDenThntiANQYYQaA1+QU\nlGjOl5u1bF1mlfbQ4AAlDe+sqwa0k7+fw6TqAFgVYQZAvSuvcOnjb7fro2+2VXmHkp/Dpmsv66CR\niZ3VtAlvsQZwdggzAOrV2p8P6rXFP500i3V8j1Yad003tWoRYlJlABoKwgyAepFTUKI3P/tJKRv3\nV2m/ICZMd17XXT06tDCpMgANDWEGQJ1yudz69Lud+uA/P1cZuTckyE9jf9tVV1zajrFiANQpwgyA\nOrMr67D+vmC9dmQertI+tG8b3XZNV4WHBplUGYCGjDAD4JyVV7i0YGm6Pv5mm1xuw9N+fstQ3Xfj\nxerOJSUA9YgwA+CcbNmVr5cWrldmdqGnzd/Prpt/00U3XN5RfowXA6CemRJmZs2apQ0bNshms2nS\npEnq0aOHZ928efP0+eefy+FwqHv37po4caIZJQI4g7IKl977cos++36HjOMnY9Ttgub648heiols\nal5xABoVr4eZ1atXa8+ePZo/f7527NihyZMna/78+ZKkwsJCvfXWW0pOTpbNZtMdd9yhjRs36uKL\nL/Z2mQCqsT3zkJ5/f60yDh4/G9Mk0KFbf9tNV8W3Yx4lAF7l9TCTmpqqxMRESVKHDh105MgRFRUV\nKSQkRAEBAQoICFBhYaGaNGmi0tJShYWFebtEAKfhcrn10bfb9MGSrVXujel9YZQeuKmnosKDTawO\nQGPl9TCTm5ur7t27e5bDw8OVm5vrCTMPPPCAEhMTFRQUpKuvvlpt27b1dokATiErt1DPv79OW/cU\neNqCAhy643fddcWlbZkQEoBpTL8B2DjhYnthYaFee+01/ec//1FISIjGjh2rrVu3qkuXLiZWCDRu\nhmEoefVe/XPxTyo7YdyYC9uG65HRvdW6BffGADCX18NMVFSUcnNzPcvZ2dmKjIyUJO3cuVNt2rTx\nXFrq27ev0tLSCDOASYpLK/TKRxu1fP3xiSEddptGX3Ghfj+kIzNbA/AJXj8SJSQkaMmSJZKktLQ0\nRUdHKzi48jp7TEyMdu7cqfLycknSpk2buMwEmGRbRoH+9PzyKkGmTXRTPfvwII1M7EyQAeAzvH5m\nJi4uTt26dVNSUpIcDoemTp2qxYsXKzQ0VImJibrjjjs0ZswY+fn5KS4uTn369PF2iUCjZhiGPvt+\np975Ik1O1/HLwMP7na+7r++hoEDTr04DQBU248SbViwoMzNTw4YNU3JysmJjY80uB7C0wuJyvfDB\neq3afMDT1iTQTw/c1FODe/P3BcA38S8WAEnSzn2HNevdVTqQV+xp6xgbpsfG9OUmXwA+jTADQN+s\n2av/9+EGlTvdnrbfDbpAt/22q/z9HCZWBgBnRpgBGrEKp0tvfLJJX6Xu9rQ1CXTo4VG9ldCztWl1\nAUBtEGaARiqnoESz56xS+t5DnrY20U018dZ+ahMdamJlAFA7hBmgEfp5d75mvrNKh46WedoG9myt\nh0bFqQnvVgJgMRy1gEZm6aq9+n8fbZDTVXl/jN1u0+3XdtPvLruAKQkAWBJhBmgkXC633vn3Zn2y\nfIenLTQ4QBNvvUQ9OrYwsTIAODeEGaARKCyp0DPvrdG6rdmetrYtQ/Xk7f3VsnmIiZUBwLkjzAAN\n3P7cIj315krtyyn0tPXv1lLjR/dWcJC/iZUBQN0gzAAN2M978vV/b/2gI0XlnrZRiZ01+ooLZbdz\nfwyAhoEwAzRQKzZm6fl5az0D4QX42fWnpN66LC7G5MoAoG4RZoAGxjAMffrdDr39eZqOzbx2XkiA\nptzeXxe2izC3OACoB4QZoAFxudx649NN+veKXZ621i1CNO2uS5lfCUCDRZgBGojSMqf+NneNVm8+\n6Gm7qF2Enry9v84LCTCxMgCoX4QZoAE4dLRMT72Zqu2Zhz1tl/WK0Z+S4hTgz0SRABo2wgxgcQfz\nizXltRTtzy3ytP1+SEeNvbor71gC0CgQZgAL27P/iKa+nqL8I5VzLNlt0r2/76mr4tuZWxgAeBFh\nBrCoLbvy9dRbK1VUUiFJ8vez67E/9FF8j9YmVwYA3kWYASxozZaDmvXuapVXuCRJTQL99OTt/XRx\nx0iTKwMA7yPMABazbG2GXpy/Xi535SAyYU0DNP2ueHWMbWZyZQBgDsIMYCGff79Tr3/yk2c5KryJ\nZtwzQDGRjCEDoPEizAAWsWDpVs396mfP8vktQzXj7ng1D2tiYlUAYD7CDODjDMPQ3K9/1sKl6Z62\nC9uGa+qdlyo0mMHwAIAwA/gwwzD09udp+mT5Dk9bz04t9OS4/goK5M8XACTCDOCz3G5D/1y8UV+l\n7Pa09b0oWhNvvYRRfQHgBIQZwAe53IZeXvijlq7e62mL79FKj/2hr/z97CZWBgC+hzAD+Biny60X\n3l+n737c52kbHBerR26Ok8NBkAGAXyPMAD6kwunWM3PXKPWn/Z624f3O1wMjesnBPEsAcEqEGcBH\nOF1u/e291Vq56YCn7eoB7XTPDRczYSQAVIMwA/iAyiCzpkqQuX5wB91+bTfZbAQZAKgOYQYwmdN1\n8qWlGy7vqHHXdCXIAEANcDchYCKXy61n561VysbjQeb6wR0IMgBQC4QZwCQul1vPvb9OKzZkedp+\nN+gCLi0BQC0RZgATuFxuPf/+On1/wtuvr73sAt35u+4EGQCoJcIM4GUul1svfLC+yjgy1yS0113X\nEWQA4GwQZgAvcrkNvTh/vZavz/S0XT2gne6+oQdBBgDOEmEG8BK329A/Fq7XsnXHg8xVA9rp3hsv\nJsgAwDkgzABeYBiGXv/kJyWvzvC0XRnfTvfeQJABgHNlyjgzs2bN0oYNG2Sz2TRp0iT16NHDs+7A\ngQMaP368nE6nunbtqunTp5tRIlBnDMPQu//erH+v2OVpS7zkfN13IyP7AkBd8PqZmdWrV2vPnj2a\nP3++nn76ac2cObPK+tmzZ+uOO+7QwoUL5XA4dODAgdPsCbCGBUvT9fG32z3Lg3rF6MGRvQgyAFBH\nvB5mUlNTlZiYKEnq0KGDjhw5oqKiIkmV/8GuXbtWQ4cOlSRNmTJFLVu29HaJQJ35ZPl2zfv6Z89y\n/24t9cjo3kwaCQB1yOthJjc3VxEREZ7l8PBw5ebmSpLy8/MVHBysmTNnavTo0Xr++ee9XR5QZ75K\n3a23PkvzLPfqHKnHx/SVn4Nb1QCgLpl+VDUMo8rn2dnZuu222zR37lxt3rxZy5cvN7E64Ox8syZD\nr368wbPc7YLmmjyunwL8HSZWBQANk9fDTFRUlOdMjCRlZ2crMjJSUuVZmpiYGMXGxsputys+Pl7b\nt28/3a4An7RiY5b+Pn+djuX0Tm2aaeod/RUUwLyuAFAfvB5mEhIStGTJEklSWlqaoqOjFRwcLEly\nOByKjY3V3r17Pevbt2/v7RKBs7Zmy0E9O3eN3L8EmXatztNTd8crOMjf3MIAoAHz+r+KcXFx6tat\nm5KSkuRwODR16lQtXrxYoaGhSkxM1KRJkzRhwgQZhqHOnTt7bgYGfF3azjzNemeVnK7KJBMT2VQz\n7olXaHCAyZUBQMNmM068acWCMjMzNWzYMCUnJys2NtbsctBI7co6rIn/738qKnVKkqIigvXXBwaq\nRbMmJlcGAA2f6TcAA1Z3IK9I015P9QSZZqGBevqeAQQZAPASwgxwDgqOlGrqa6kqOFomSQoO8tNT\nd8WrVYsQkysDgMaDMAOcpaKSCk17I1X78yoHfQzws2vK7f11QUyYyZUBQONyVmHG7XbXdR2ApZRV\nuPR/b/+gXVlHJEl2u02Pj+mr7h1amFwZADQ+NQozixYt0rx58+R0OnXzzTdr2LBhev/99+u7NsAn\nuVxuPfPeGqXtzPO0/XFEL/Xv3srEqgCg8apRmFmwYIFGjBihpUuXqlOnTkpOTtZXX31V37UBPscw\nDL384Qb9kHZ8AtRx13RTYr/zTawKABq3GoWZwMBABQQEaPny5brqqqtkt3OrDRqnd77YrKWr93qW\nfz+ko24c0tHEigAANU4lTz31lNatW6d+/fpp/fr1Ki8vr8+6AJ+z6NttWrTs+PQaw/udr1t/29XE\nigAAUg3DzLPPPqu2bdvq1VdflcPh0L59+/TUU0/Vd22Az1i6ao/+9cVmz/Kl3VvqgZt6ymazmVgV\nAECq4XQGb7zxhiZPnuxZvuaaa+qtIMDXpP60X/9Y+KNnuXuH5nrsD33lcHC5FQB8QY2Oxg6HQ6mp\nqSorK5Pb7fZ8AA3dT9tz9cwJE0de0DpMT47rrwB/h7mFAQA8anRm5sMPP9S7774rwzBks9k8j1u2\nbKnv+gDT7Mg8pP97+wdVOCuDe6sWIZp+96UKacIM2ADgS2oUZtauXVvfdQA+JSunUNPfWKmSssr5\nliLOC9SMu+MVHhpkcmUAgF+r0WWmw4cP669//asee+wxSdI333yj/Pz8ei0MMEve4RJNeT1Vhwor\n51sKaeKvGXcPUMvmzLcEAL6oRmHmySefVKtWrZSRkSFJKi8v1xNPPFGvhQFmKCwu17TXU5WdXyxJ\nCvB3aNodl6ptq/NMrgwAcDo1CjP5+fkaO3as/P0r7xW48sorVVpaWq+FAd5WWu7UjLd+0J4DRyVJ\nDrtNE2+9RBe1jzC5MgBAdWr83tKKigrPmBq5ubkqLi6ut6IAb3O63PrrnDXasvv45dM/JcWp70XR\nJlYFAKiJGt0AfMstt+imm25STk6O7r33Xv30009Vxp0BrMztNvT3Beu1ZstBT9td13XX5X3amFgV\nAKCmahRmrr76avXu3Vvr169XQECAZsyYoaioqPquDah3hmHorc83adnaTE/byMTO+t2gDiZWBQCo\njWrDzCeffHJSW1lZmVJSUiRJ119/ff1UBXjJh8nb9Nl3Oz3LV8a30x+uvNDEigAAtVVtmFmxYoUk\nqaCgQD///LN69uwpl8uljRs3Ki4ujjADS/s6dbfe++r4wI8JF7fWvTdezHxLAGAx1YaZZ555RpL0\n0EMPaenSpQoKqhwwrLCwUE8++WT9VwfUkxUbs/Tqxxs8yz07tdCjt/SWw06QAQCrqdG7mbKysjxB\nRpKaNm2qrKyseisKqE8b0nP07Ny1nvmWOrZppkm39ZO/H/MtAYAV1egG4E6dOikpKUlxcXGy2+3a\nsGGDzj///PquDahz2zIKNPOdH+R0Vc63FBPZVNPvvFTBQcy3BABWVaMw85e//EUpKSlKT0+XYRi6\n6667dNlll9V3bUCdyjh49Jf5llySpOZhQZpxT7zCmgaaXBkA4FxUe5lp8+bNkqSVK1fKbrfrwgsv\n1EUXXaSAgAD98MMPXikQqAu5h0o09fVUHSkqlySFBvtrxt3xigoPNrkyAMC5qvbMzKeffqquXbvq\nlVdeOWmdzWZTfHx8vRUG1JUjReWa+nqKcg+VSJICAxyaeuelOr8l8y0BQENgMwzDMLuIc5GZmalh\nw4YpOTlZsbGxZpcDH1NS5tSUf6Zo694CSZKfw6Ypt1+q3hcy6CMANBQ1ejfTqlWrdOONN6pnz57q\n1auXRo0apR9//LG+awPOSYXTpb+8s8oTZGw26ZGbexNkAKCBqfENwE888YT69OkjwzC0Zs0aTZ8+\n/ZQjBAO+wOU29Nz76/Rjeo6n7d4bL9agOM7eAUBDU6MzM82aNVN8fLwCAgIUGBiohIQERUczmzB8\nk2EYevXjDVqx4fhYSH+48kJdPaC9iVUBAOpLjc7M9OzZU++8844GDhwot9utlStXqkOHDsrIyJAk\ntWnD7MLwHe99tUVLVu7xLP/usgs0MrGziRUBAOpTjcLM559/Lkl67733PG2GYejrr7+WzWZTcnJy\n/VQH1NLiZdv1YfI2z/KQPrG643fdmW8JABqwasNMYWGhPvroI33zzTeSpA8++EAffPCB2rVrp6lT\np6pFixZeKRKoiaWr9ujtz9M8y/26ttRDo+JkZ74lAGjQqr1nZurUqcrLy5Mk7dq1Sy+88IImTJig\nhIQEzZw50ysFAjWR+tN+/WPh8XfYdbuguR4f21d+jhrdFgYAsLBqz8xkZGTo+eeflyQtWbJEV155\npQYMGCBJ+uKLL+q/OqAGNm7P0d/eW+OZOPKC1mGacnt/BfozcSQANAbV/tsaHHx8qPdVq1bp0ksv\n9SxzDwJ8wbaMAj399vGJI1u1CNH0uy9VSBMmjgSAxqLaMONyuZSXl6e9e/dq/fr1SkhIkCQVFRWp\npKTEKwUCp3OqiSP/754BCg8NMrkyAIA3VXuZ6a677tLVV1+t0tJSPfjggwoLC1NpaalGjx6tkSNH\neqtG4CQ5BaeeODI6gokjAaCxqTbMDB48WP/73/9UVlampk2bSpKCgoL02GOPaeDAgWf9pLNmzdKG\nDRtks9k0adIk9ejR46RtnnvuOf34449V3g4OSFLB0VJNeW2FZ+LIoACHpjFxJAA0WmccZ8bf31/+\n/lXvPzjp3Xb5AAAa10lEQVSXILN69Wrt2bNH8+fP144dOzR58mTNnz+/yjY7duzQmjVrTnpe4Ghx\nuaa+lqp9OUWSKieOnHRbP3VpG2FyZQAAs3j9faupqalKTEyUJHXo0EFHjhxRUVFRlW1mz56t8ePH\ne7s0+Lji0gpNez1Vu/cfkSTZ7TY9Pqav4rowcSQANGZeDzO5ubmKiDj+X3R4eLhyc3M9y4sXL1b/\n/v3VunVrb5cGH1Za7tSMt37QtoxDkipnwP5TUpzie/B7AgCNnekjihmG4fn88OHDWrRokcaNGyfD\nMKqsQ+NV4XRp1jurlbYzz9N2/+97akgf5gQDAJgQZqKioqqcicnOzlZkZKQkaeXKlSooKNAtt9yi\nP/7xj9qyZYtmz57t7RLhQ1wut56Zu1brtmZ72u74XTddGd/OvKIAAD7F62EmISFBS5YskSSlpaUp\nOjraMzjfFVdcoS+++ELz58/Xyy+/rK5du2rChAneLhE+wu029OKC9Ur9ab+nbfRvuuj6wR1NrAoA\n4GtqNGt2XYqLi1O3bt2UlJQkh8OhqVOnavHixQoNDfXcGAwYhqFXF23UsrWZnrYbLu+opN90MbEq\nAIAvshkWvzElMzNTw4YNU3JysmJjY80uB3XAMAy9/XmaPlm+w9N2VXw73ff7i5lGAwBwEtNvAAZO\nZBiG3v335ipBZkifWN17I0EGAHBqhBn4DMMw9N5XW/Txt9s9bfE9WunhUXGy2wkyAIBTI8zAZ7y/\nZKs+TN7mWe7fraUe+0NfORz8mgIATo9XCfiED/6zVfP/u9WzfEnXaD0x9hL5+/ErCgCoHq8UMN3C\npel6f8nPnuW+F0Vr4q0EGQBAzfBqAVN99M02vffVFs9y7y5RvwQZh4lVAQCshDAD0yxetl3v/nuz\nZ7lXp0hNGtdPAf4EGQBAzRFmYIpPlu/Q25+neZYv7thCk2/vp0CCDACglrw+AjDwyfLteuuz40Gm\ne4fmmnJ7fwUF8OsIAKg9Xj3gVQuXple5R6bbBc019Y5LFRTIryIA4OzwCgKvMAxD85b8rAX/Tfe0\nVQaZ/mpCkAEAnANeRVDvjk1RcOLIvr06RWry7f24tAQAOGe8kqBeGYahNz7dpM+/3+lpOzaODO9a\nAgDUBcIM6o3bbejVRRv1depuT9ul3Vvq8TF9GUcGAFBnCDOoFy63oX8sXK/k1RmetoE9W+vRW/rI\nj7mWAAB1iDCDOlfhdOmZuWuV+tN+T9uQPrF6eFQck0YCAOocYQZ1qri0Qn95Z5U2bMv1tA3vd74e\nGNFLDrvNxMoAAA0VYQZ15khRuZ56M1Xpew952q4b1EG3X9tNdoIMAKCeEGZQJ3IPlWjq6ynKOFjo\naRtz1UUaMayTbDaCDACg/hBmcM725RRqymspyikokSTZbNJ9N16sqwa0N7kyAEBjQJjBOdmeeUjT\n30jV4cJySZLDbtOjo/vosrgYkysDADQWhBmctQ3pOfrLu6tUXOqUJAUGODTx1kvU58JokysDADQm\nhBmclW/W7NVLC36Uy21IkkKa+GvaHZfqovYRJlcGAGhsCDOoFcMwtHBpuuZ+/bOnrXlYkKbfFa92\nrc4zsTIAQGNFmEGNOV1uvfLRBv131V5PW7tW52nanZeqRbMmJlYGAGjMCDOokeLSCv11zhqt25rt\naevZqYUm3tpPIU38TawMANDYEWZwRnmHSzTjzR+0M+uwp21o3zZ6cEQv+fsxPQEAwFyEGVRrV9Zh\nzXjrB+UeKvG0JQ3votFXdGEwPACATyDM4LRSf8rSc++vU1m5S5Jkt9v0wE099Zv+bU2uDACA4wgz\nOIlhGPoweZve+2qLp61JoJ+eGNuXMWQAAD6HMIMqyipcemnBen23fp+nrWXzYE25vb/Ob8lbrwEA\nvocwA4/8I6Wa+a8fqsx63b1Dc00Ye4nCmgaaWBkAAKdHmIEkacuufM2es0r5R8o8bVdc2lb33HAx\n71gCAPg0wkwjZxiGvlyxS298uskzNYHdJt15XQ9dM7A971gCAPg8wkwjVlbh0isfbdA3azI8baHB\nAXrsD30U1yXKxMoAAKg5wkwjdSCvSLPeWV1lILwOsWGadGs/RUUEm1gZAAC1Q5hphFZtPqAX3l+n\nwpIKT9uwS9rovt/3VKC/w8TKAACoPcJMI1LhdGvOl5v1yfIdnjY/h013X99DV8a34/4YAIAlEWYa\niQN5RXpm7poqb7tuHhakCbdeogvbRphYGQAA58aUMDNr1ixt2LBBNptNkyZNUo8ePTzrVq5cqRde\neEEOh0Pt27fXzJkzzSixQUnZmKWXFqxXUanT09b3omj9KSmO8WMAAJbn9TCzevVq7dmzR/Pnz9eO\nHTs0efJkzZ8/37N+2rRpeu+99xQVFaWHH35Y3333nQYNGuTtMhuE0nKn/vV5mr5M2e1pc9htuvW3\nXXXdoA6y27msBACwPq+HmdTUVCUmJkqSOnTooCNHjqioqEghISGSpEWLFnk+j4iI0KFDh067L5ze\n9sxDem7eWmVmF3raosKb6PExfdWFy0oAgAbE60O75ubmKiLi+ItpeHi4cnNzPcvHgkx2drZSUlI0\nePBgb5doaS63oQ+T0/Xnv39XJcjE92ilv4+/nCADAGhwTL8B2DCMk9ry8vJ03333afr06QoLCzOh\nKms6kFekFz5Yp8278j1tQQEO3XldD/2m//m8WwkA0CB5PcxERUVVOROTnZ2tyMhIz3JhYaHuuusu\nPfroo4qPj/d2eZZkGIb+88NevfXZJpWUHb/Jt0vbcI0f3VutWzQ1sToAAOqX1y8zJSQkaMmSJZKk\ntLQ0RUdHKzj4+Iizs2fP1rhx45SQkODt0izpYH6xpr6Wqpc//NETZOx2m0ZfcaH++sBAggwAoMHz\n+pmZuLg4devWTUlJSXI4HJo6daoWL16s0NBQDRw4UJ999pn27t2rhQsXymaz6dprr9WIESO8XabP\nc7sNfZWyS+/8e7NKy12e9pjIEI0f3Uedzw83sToAALzHZpzqphULyczM1LBhw5ScnKzY2Fizy/GK\nrNxCvbTgR6XtzPO02W3SdYM76pYrL2RKAgBAo2L6DcCouQqnS4u+3a6FS9NV7nR72ttEh+rhUb14\npxIAoFEizFjEhm05evXjjdqXc/zt1na7TSOGdtKo4Z3l78fZGABA40SY8XEFR0v19mdpWrYus0p7\nx9gwPTCilzrGNjOpMgAAfANhxof9vDtf099IrTKnUnCQn8ZcdZGuGtBeDqYjAACAMOPLFixNrxJk\nBvWK0R3XdVfEeUEmVgUAgG8hzPiw+B6ttG5rtmIiQ3TXdT0U1yXK7JIAAPA5hBkf9pv+bXV571gF\n8FZrAABOy+sjAKN2CDIAAFSPMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMA\nACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyN\nMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMA\nACyNMAMAACyNMAMAACyNMAMAACyNMAMAACzNz4wnnTVrljZs2CCbzaZJkyapR48ennUpKSl64YUX\n5HA4NGjQIN1///1mlAgAACzC62dmVq9erT179mj+/Pl6+umnNXPmzCrrZ86cqZdfflkffPCBVqxY\noR07dni7RAAAYCFeDzOpqalKTEyUJHXo0EFHjhxRUVGRJCkjI0PNmjVTdHS0bDabBg8erJUrV3q7\nRAAAYCFeDzO5ubmKiIjwLIeHhys3N/eU6yIiIpSdne3tEuGjtm3bpjZt2igkJERt2rTx6lm7vLw8\njRo1SnFxcWrTpo169+6tUaNGKT8//7Tb9uvXT9ddd51uuOEG9evX76TtT9zudPvyptr0EQB8iSn3\nzJzIMIyzWofGZ+jQocrMzJQkFRcX6/LLL1dGRoZXnvv+++/XwoULPcuZmZlav369JGnBggXVbnvM\n6tWrq2x/4na/XmeG2vQRAHyJ18/MREVFec7ESFJ2drYiIyM963JycjzrDh48qKioKG+XCB/16zME\n3jxjsGvXrhq3n27bX6/79XbVfZ031KaPAOBLvB5mEhIStGTJEklSWlqaoqOjFRwcLEmKiYlRUVGR\nsrKy5HQ6tWzZMg0cONDbJcJHnXgJ8lTL9al9+/Y1bj/dtr9e9+vtqvs6b6hNHwHAl3j9MlNcXJy6\ndeumpKQkORwOTZ06VYsXL1ZoaKgSExM1bdo0jR8/XpJ0zTXXqG3btt4uET5q2bJluvzyy5Wfn6+I\niAgtW7bMa8/96quvSpLS09OVm5uryMhIderUydN+qm137dql1q1by2azad++fWrfvn2V7U/c7tfr\nzFCbPgKAL7EZFr8xJTMzU8OGDVNycrJiY2PNLgcAAHgZIwADAABLI8wAAABLI8wAAABLI8wAAABL\nI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wA\nAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABL\nI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wA\nAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABLI8wAAABL8/P2EzqdTk2YMEFZWVlyOBya\nNWuWYmNjq2zz5Zdf6l//+pccDof69++vRx55xNtlAgAAi/D6mZkvvvhCYWFhev/993Xvvffqueee\nq7K+tLRUzz33nObMmaP58+crNTVVO3bs8HaZAADAIrweZlJTU5WYmChJGjBggNatW1dlfVBQkD7/\n/HM1adJEktSsWTMdOnTI22UCAACL8HqYyc3NVUREhCTJZrPJbrfL6XRW2SY4OFiStHXrVmVlZalX\nr17eLhMAAFhEvd4z8+GHH+qjjz6SzWaTJBmGoY0bN1bZxu12n/Jrd+/erT//+c967rnn5HA4Tvsc\nLpdLknTgwIE6qhoAAJipZcuW8vOreUSp1zAzYsQIjRgxokrbxIkTlZubqy5dunjOyPy64AMHDuiP\nf/yjnnnmGXXp0qXa58jJyZEk3XLLLXVYOQAAMEtycvJJbw6qjtffzZSQkKCvv/5aCQkJ+uabb9S/\nf/+Ttpk8ebKmTZumCy+88Iz76969u+bNm6fIyMhqz+AAAABraNmyZa22txmGYdRTLafkdrs1efJk\n7dmzR4GBgZo9e7aio6P1+uuvq3///goLC9MNN9ygHj16yDAM2Ww2jRs3TkOGDPFmmQAAwCK8HmYA\nAADqEiMAAwAASyPMAAAASyPMAAAAS7NkmHE6nfrzn/+s0aNHa8yYMcrMzDxpmy+//FIjRoxQUlKS\nXnjhBROqrL1Zs2YpKSlJN998s3766acq61JSUjz9eeWVV0yq8NxU17+VK1dq1KhRGj16tCZPnmxS\nhWevur4d89xzz2nMmDFerqxuVNe/AwcOaPTo0Ro5cqSmT59uToHnqLr+zZs3T0lJSbrllls0a9Ys\nkyo8e+np6Ro+fLjmzZt30rqGcFyprn9WP65I1ffvGKseW6rrW62PK4YFLV682JgxY4ZhGIbxv//9\nz/jTn/5UZX1JSYkxdOhQo7i42DAMwxgxYoSxfft2r9dZG6tWrTLuuecewzAMY/v27caoUaOqrL/6\n6quNAwcOGG632xg9erTP9+fXztS/3/zmN8bBgwcNwzCMhx56yFi+fLnXazxbZ+rbsfakpCRjzJgx\n3i7vnJ2pfw8//LCxdOlSwzAMY8aMGcb+/fu9XuO5qK5/R48eNYYMGWK43W7DMAzj9ttvNzZs2GBK\nnWejuLjYGDNmjDFlyhRj7ty5J623+nHlTP2z8nHFMM7cP8Ow7rHlTH2r7XHFkmdmGuL8Tif2qUOH\nDjpy5IiKiookSRkZGWrWrJmio6Nls9k0ePBgrVy50sxya626/knSokWLFBUVJUmKiIjw+Z/Xic7U\nN0maPXu2xo8fb0Z556y6/hmGobVr12ro0KGSpClTptR6fAizVde/gIAABQQEqLCwUE6nU6WlpQoL\nCzOz3FoJDAzUm2++6fnbOlFDOK5U1z/J2scV6cz9k6x7bKmub2dzXLFkmGmI8zud2CdJCg8PV25u\n7inXRUREKDs72+s1novq+idJISEhkqTs7GylpKRo8ODBXq/xbJ2pb4sXL1b//v3VunVrM8o7Z9X1\nLz8/X8HBwZo5c6ZGjx6t559/3qwyz1p1/QsICNADDzygxMREDRs2TBdffLHatm1rVqm1ZrfbFRAQ\ncMp1DeG4Ul3/JGsfV6Qz98/Kx5bq+nY2xxWvjwBcW96Y38kXGdUM/1PdOqs4VR/y8vJ03333afr0\n6Zb67/fXTuzb4cOHtWjRIr3zzjvav39/g/vZGYah7Oxs3XbbbWrdurXuvvtuLV++3HIvGic6sX+F\nhYV67bXX9J///EchISEaO3astm7desZpVqyoIfxunkpDOa78WkM8thxzNscVnw8z3pjfyRdERUVV\n+W8+OztbkZGRnnXH5qCSpIMHD1Z72tEXVdc/qfJF46677tKjjz6q+Ph4M0o8a9X1beXKlSooKNAt\nt9yisrIyZWRkaPbs2ZowYYJZ5dZadf0LDw9XTEyMZw6V+Ph4bd++3VJhprr+7dy5U23atPG8CPbt\n21dpaWmWOKacSUM4rpyJlY8rZ9IQji2nczbHFUteZjo2v5OkOpnfyRckJCRoyZIlkqS0tDRFR0d7\nLpXFxMSoqKhIWVlZcjqdWrZsmQYOHGhmubVWXf+kyuu+48aNU0JCglklnrXq+nbFFVfoiy++0Pz5\n8/Xyyy+ra9euljvYVNc/h8Oh2NhY7d2717O+ffv2ptV6Ns70t7dz506Vl5dLkjZt2mSpy0zVaQjH\nlTOx8nHlTBrCseV0zua4YsnpDBrq/E7PP/+8Vq1aJYfDoalTp2rz5s0KDQ1VYmKi1qxZo2effVaS\ndOWVV+q2224zt9izcLr+DRw4UP369VOvXr08P69rr732pDNyvqy6n90x+/bt08SJEzVnzhwTKz07\n1fVv7969mjBhggzDUOfOnfXUU0+ZXW6tVde/hQsX6uOPP5afn5/i4uL05z//2exyaywtLU2zZ89W\nVlaW/Pz8FB0draFDhyo2NrZBHFeq619DOK6c6ed3jBWPLWfqW22PK5YMMwAAAMdY8jITAADAMYQZ\nAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAFXk5uZqwoQJuuGGG3TLLbfohhtuOKvxK2bN\nmqXNmzefUy2vv/66li9ffk77kCpH9U1KStLhw4fPeV/HPPTQQ1qxYkWd7Q/A2WOcGQBVjBw5Ujfd\ndJNGjhwpqXJum9tuu00PPfSQhg8fbnJ1Z+eee+7RrbfeqgEDBtTZPg8dOqSRI0fq008/VZMmTeps\nvwBqz+fnZgLgPampqfLz8/MEGUlq3ry5Fi9e7Jn/bOLEifL399fu3bv17LPPauPGjXrzzTcVGBgo\nl8ulv/3tb2rdurXGjBmj+++/Xw6HQ6+//rpatmyp7du3y9/f37P9McXFxXr00Ud19OhROZ1ODRky\nRPfcc48mTpyoPn36KCIiQv/6179ks9lUXFysbdu2acOGDTpy5IimTZumgoICHT16VOPGjdM111xT\npU9btmzR/v37PUFm6NChuvnmm/X9998rJydHTzzxhObPn6+dO3fq/vvv1/XXX6+JEyeqWbNm2rlz\np7Zv365HH31U33zzjbZu3ao+ffpo+vTpatasmYYMGaIPP/xQY8eO9cJPB8DpcJkJgMe2bdvUvXv3\nk9p/PZFraWmp5syZo6ioKB09elQvvvii3n33XQ0aNEhz58496es3bNigRx99VPPnz5fNZtP3339f\nZX1KSopcLpfmzp2rDz74QMHBwVVmAR46dKjee+89zZkzR7GxsZoyZYok6cUXX9SgQYP0zjvvaO7c\nuXrppZdUUFBQZd/ff/+9LrvssiptzZs315w5c9SrVy/NmTNH//znP/X000/r3Xff9WyTn5+v1157\nTQ8++KBmzJih6dOn66OPPtLixYtVWFgoSRowYMBJfQHgfZyZAeDhcDjkcrk8ywsXLtTnn3+u8vJy\ntWrVSi+++KIkKS4uzrNNRESEHn/8cRmGodzcXPXq1euk/Xbo0EHh4eGSKic4/PW9K71799ZLL72k\nRx55RIMGDdKIESNks9lO2s9bb72l0NBQ3XTTTZKkH374QZs2bdKiRYskSQEBAcrMzPQ8lyTt379f\nHTp0qLKfY/VHR0erZcuWkqSWLVvq6NGjVWo61t6hQwc1bdpUUuWMvkePHlXTpk3VunVr7du37/Tf\nUABeQZgB4NGlSxd9/PHHnuWRI0dq5MiRWrVqlf7+97972v39/SVJTqdTjzzyiD799FO1adNG8+bN\n06ZNm07ar8PhqLL861v1IiIi9Nlnn2n9+vVaunSpbrzxRn3yySdVtlm5cqX++9//VjnzExAQoGnT\npqlbt2616ueJZ5pOrO3Euk5s//U23GoI+BYuMwHw6Nu3r8LDw/X666972ioqKrRixQoFBQWdtH1R\nUZEcDodat26tsrIyJScnq7y8vNbPu2LFCn377beKi4vTY489ppCQEOXl5XnWHzx4UE8//bRefPHF\nKkGkT58++vLLLyVVXvp66qmn5Ha7q+y7VatW2r9/f61rqomsrCzFxMTUy74B1BxhBkAVr776qvLy\n8nT99ddrzJgxSkpKUmlpqZ599tmTtg0LC9M111yj3//+9xo/frzuvPNO/fDDD1qyZMkpLxOdTvv2\n7fX222/rD3/4g8aOHauBAweqVatWnvWvvPKKioqK9Pjjj2vMmDEaO3asdu/erQcffFB79uzR6NGj\nNWbMGHXt2lV2e9XD2mWXXVblvpba1HUqJ359SkrKSffjAPA+3poNoMG79957NWbMGCUkJNTZPgsK\nCpSUlKTFixcrODi4zvYLoPYIMwAavJycHD300EP65z//qbCwsDrZ50MPPaSbb75Z8fHxdbI/AGeP\nMAMAACyNe2YAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAICl/X8Fuim4qqHHmAAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8812f74fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_log_reg(x='Grain size (mm)', y='Spiders', data=df, clf=clf, xmin=0, xmax=1.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis testing\n",
    "\n",
    "To test if *Grain size* is a significant factor, we use the [**likelihood ratio test**](https://en.wikipedia.org/wiki/Logistic_regression#Evaluating_goodness_of_fit).\n",
    "\n",
    "We calculate the likelihood of the model with the grain size (the alternative model):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def log_reg_null_model(y):\n",
    "    clf = sklearn.linear_model.LogisticRegression(C=1e12)\n",
    "    clf.fit(np.zeros_like(y).reshape(-1, 1), y)\n",
    "    return clf\n",
    "\n",
    "clf0 = log_reg_null_model(df['Spiders'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The likelihood ratio test operates by calculating the test statistic $D$ from the likelihoods of the null and alternative models:\n",
    "$$\n",
    "D = -2 \\log{ \\frac{L(H_0)}{L(H_1)} }\n",
    "$$\n",
    "The test statistic is then approximately chisquare distributed.\n",
    "\n",
    "*scikit-learn* has a log-loss function that can help us do that. \n",
    "The log-loss is defined as the negative log-likelihood, so we can rewrite:\n",
    "$$\n",
    "D = 2 (-\\log{L(H_0)} + \\log{L(H_1)}) \\Rightarrow \\\\\n",
    "D = 2 (logloss(H_0) - logloss(H_1))\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sklearn.metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def log_reg_lik_ratio_test(X, Y, clf0, clf1, df=1):\n",
    "    if X.ndim == 1:\n",
    "        X = X.reshape(-1, 1)\n",
    "    y_prob0 = clf0.predict_proba(X)\n",
    "    loss0 = sklearn.metrics.log_loss(Y, y_prob0, normalize=False)\n",
    "    y_prob1 = clf1.predict_proba(X)\n",
    "    loss1 = sklearn.metrics.log_loss(Y, y_prob1, normalize=False)\n",
    "    D = 2 * (loss0 - loss1)\n",
    "    return scipy.stats.distributions.chi2.sf(D, df=df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.033243767135570736"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_reg_lik_ratio_test(df['Grain size (mm)'], df['Spiders'].astype(np.float64), clf0, clf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "John indeed reports 0.033.\n",
    "\n",
    "Note that the log-loss calculation in equivalent to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15.3157739197 15.3157739197\n"
     ]
    }
   ],
   "source": [
    "_ = clf.predict_proba(df['Grain size (mm)'].reshape(-1, 1))\n",
    "df['prob_absent'], df['prob_present'] = _[:,0], _[:,1]\n",
    "lik = df.loc[df['Spiders'], 'prob_present'].prod() * df.loc[~df['Spiders'], 'prob_absent'].prod()\n",
    "print(\n",
    "    -np.log(lik), \n",
    "    sklearn.metrics.log_loss(\n",
    "        df['Spiders'], \n",
    "        clf.predict_proba(df['Grain size (mm)'].reshape(-1, 1)), \n",
    "        normalize=False\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Second example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Location</th>\n",
       "      <th>Latitude</th>\n",
       "      <th>Mpi90</th>\n",
       "      <th>Mpi100</th>\n",
       "      <th>p, Mpi100</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Port Townsend, WA</td>\n",
       "      <td>48.1</td>\n",
       "      <td>47</td>\n",
       "      <td>139</td>\n",
       "      <td>0.748</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Neskowin, OR</td>\n",
       "      <td>45.2</td>\n",
       "      <td>177</td>\n",
       "      <td>241</td>\n",
       "      <td>0.577</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Siuslaw R., OR</td>\n",
       "      <td>44.0</td>\n",
       "      <td>1087</td>\n",
       "      <td>1183</td>\n",
       "      <td>0.521</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Umpqua R., OR</td>\n",
       "      <td>43.7</td>\n",
       "      <td>187</td>\n",
       "      <td>175</td>\n",
       "      <td>0.483</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Coos Bay, OR</td>\n",
       "      <td>43.5</td>\n",
       "      <td>397</td>\n",
       "      <td>671</td>\n",
       "      <td>0.628</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            Location  Latitude  Mpi90  Mpi100  p, Mpi100\n",
       "0  Port Townsend, WA      48.1     47     139      0.748\n",
       "1       Neskowin, OR      45.2    177     241      0.577\n",
       "2     Siuslaw R., OR      44.0   1087    1183      0.521\n",
       "3      Umpqua R., OR      43.7    187     175      0.483\n",
       "4       Coos Bay, OR      43.5    397     671      0.628"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = \"\"\"Location\tLatitude\tMpi90\tMpi100\tp, Mpi100\n",
    "Port Townsend, WA\t48.1\t47\t139\t0.748\n",
    "Neskowin, OR\t45.2\t177\t241\t0.577\n",
    "Siuslaw R., OR\t44\t1087\t1183\t0.521\n",
    "Umpqua R., OR\t43.7\t187\t175\t0.483\n",
    "Coos Bay, OR\t43.5\t397\t671\t0.628\n",
    "San Francisco, CA\t37.8\t40\t14\t0.259\n",
    "Carmel, CA\t36.6\t39\t17\t0.304\n",
    "Santa Barbara, CA\t34.3\t30\t0\t0\n",
    "\"\"\"\n",
    "df = pd.read_table(io.StringIO(data))\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAFhCAYAAACCiIhHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHrBJREFUeJzt3X9Y1fX9//HHG1AK4xIwOCSol5E/gjQxWzOdpLC1HNVV\nV+4iFdRRW5muRT+muNRP6qSszGnt2q71Y6ZJadLSVWydpq6JSq5wUi7DqaABHjNREhE53z+6PF9/\ncDxAvM+BF/fbdXldvl7v4znP07u3D1/v94vXy3K73W4BAAAjBAW6AAAA0HYIdgAADEKwAwBgEIId\nAACDEOwAABikwwR7Q0ODKioq1NDQEOhSAABotzpMsFdWVio1NVWVlZWBLgUAgHarwwQ7AADwjWAH\nAMAgBDsAAAYh2AEAMAjBDgCAQQh2AAAMQrADAGAQgh0AAIMQ7AAAGIRgBwDAIAQ7AAAGIdgBADAI\nwQ4AgEEIdgAADEKwAwBgEIIdAACDEOwAABiEYAcAwCAEOwAABiHYAQAwiO3B/vnnn+uHP/yhVq5c\necGxzZs3a9y4ccrIyNALL7xgdykAABjP1mA/ceKE5s+fr+HDhzd5fMGCBVq2bJlWrVqlf/3rXyor\nK7OzHAAAjGdrsIeGhupPf/qTYmJiLjhWXl6uiIgIORwOWZallJQUbdmyxc5yAAAwnq3BHhQUpK5d\nuzZ5zOVyKSoqytOOiopSdXW1neUAAGC8djN5zu12B7oEAAA6vIAFe0xMjA4dOuRpV1VVNXnLHgAA\nNF/Agj0uLk61tbU6ePCgGhoatGHDBo0cOTJQ5QAAYIQQO9+8tLRUeXl5OnjwoEJCQlRYWKgxY8Yo\nPj5eaWlpmjNnjnJyciRJ6enp6tOnj53lAABgPMvdQR5uV1RUKDU1VU6nU/Hx8YEuBwCAdqndTJ4D\nAADfHcEOAIBBCHYAAAxCsAMAYBCCHQAAgxDsAAAYhGAHAMAgBDsAAAYh2AEAMAjBDgCAQQh2AAAM\nQrADAGAQgh0AAIMQ7AAAGIRgBwDAIAQ7AAAGIdgBADAIwQ4AgEEIdgAADEKwAwBgEIIdAACDEOwA\nABiEYAcAwCAEOwAABiHYAQAwCMEOAIBBCHYAAAxCsAMAYBCCHQAAgxDsAAAYhGAHAMAgBDsAAAYh\n2AEAMAjBDgCAQQh2AAAMQrADAGAQgh0AAIMQ7AAAGIRgBwDAIAQ7AAAGIdgBADAIwQ4AgEFC7P6A\nhQsXqqSkRJZlKTc3V4MGDfIcW7lypdatW6fg4GBdc801mjlzpt3lAABgNFuDvbi4WPv27VN+fr7K\nyso0a9Ys5efnS5KOHz+uF198UU6nU5ZlKTs7Wzt27NDgwYPtLAkAAKPZeiu+qKhIaWlpkqSEhATV\n1NSotrZWktS1a1d17dpVx48fV0NDg+rq6tS9e3c7ywEAwHi2BrvL5VJUVJSnHRkZKZfLJenbYH/g\ngQeUlpam1NRUDR48WH369LGzHAAAjOfXyXNut9vz++PHj+sPf/iD/va3v8npdKqkpET//e9//VkO\nAADGsTXYY2JiPCN0SaqurlZ0dLQkac+ePerVq5e6d++ukJAQDRs2TKWlpXaWAwCA8WwN9hEjRqiw\nsFCSVFpaKofDobCwMElSXFyc9uzZo/r6eknSzp07uRUPAMB3ZOus+OTkZCUlJSkjI0PBwcGaPXu2\nCgoKFB4errS0NGVnZyszM1MhISFKTk7WddddZ2c5AAAYz3Kf/eC7HauoqFBqaqqcTqfi4+MDXQ4A\nAO0SK88BAGAQgh0AAIMQ7AAAGIRgBwDAIAQ7AAAGIdgBADCI7du2AgDQXmz6uEKrnbu1v+qYejvC\nNS61n0Ylm/Uj1AQ7AKBT2PRxhRat2O5p7/2yxtM2Kdy5FQ8A6BRWO3e3qL+jItgBAJ3C/qpjTfaX\ne+nvqAh2AECn0NsR3mR/Ly/9HRXBDgDoFMal9mtRf0fF5DkAQKdwZoLcaudulVcdUy9mxQMA0LGN\nSo43LsjPx614AAAMQrADAGAQgh0AAIMQ7AAAGIRgBwDAIAQ7AAAGIdgBADAIwQ4AgEEIdgAADEKw\nAwBgEIIdAACDEOwAABiEYAcAwCAEOwAABiHYAQAwSLP2Y9+1a5e+/vprud1uT9/w4cNtKwoAALSO\nz2CfPn26du3apdjYWE+fZVkEOwAA7ZDPYD9w4ID+/ve/+6MWAADwHfl8xt63b1/V19f7oxYAAPAd\n+RyxBwUF6Sc/+YkGDx6s4OBgT/9TTz1la2EAAKDlfAb7jTfeqBtvvNEftQAAgO/IZ7Dfcccdqqio\n0KeffirLspSUlKSePXv6ozYAANBCPp+xr1q1SllZWfrrX/+qdevWKTMzUwUFBf6oDQAAtJDPEftf\n/vIXvfvuuwoNDZUkffPNN5oyZYruuOMO24sDgI5q08cVWu3crf1Vx9TbEa5xqf00Kjk+0GWhE/AZ\n7CEhIZ5Ql6SwsDB16dLF1qIAoCPb9HGFFq3Y7mnv/bLG0ybcYTefwR4bG6t58+Z5JtB9+OGHuuKK\nK2wvDAA6qtXO3V77CXbYzWewz5s3T6+++qrWrl0ry7J07bXXKjMz0x+1AUCHtL/qWJP95V76gbbk\nNdjdbrcsy1JoaKjuuecef9YEAB1ab0e49n5Zc0F/L0d4AKpBZ+M12CdNmqTly5crMTFRlmV5+s8E\n/meffdasD1i4cKFKSkpkWZZyc3M1aNAgz7HKykrl5OSooaFBiYmJmjt3buu/CQC0E+NS+53zjP3s\nfsBuXoN9+fLlkqStW7eqe/fu5xwrLy9v1psXFxdr3759ys/PV1lZmWbNmqX8/HzP8by8PGVnZys1\nNVXz5s1TZWXlOZvNAEBHdOY5+mrnbpVXHVMvZsXDjy76jL2xsVHTpk3T8uXLPSP1U6dOaerUqVq3\nbp3PNy8qKlJaWpokKSEhQTU1NaqtrVW3bt3kdru1fft2LV68WJL0+OOPt8HXAYD2YVRyPEGOgPAa\n7OvXr9fSpUu1b98+JSYmevZiDwoK0siRI5v15i6XS9dcc42nHRkZKZfLpW7duumrr75SWFiYFixY\noE8//VTDhg1TTk7Od/w6AAB0bl6DPT09Xenp6Vq6dKmmT5/eJh925h8HZ35fXV2tyZMnq2fPnvr5\nz3+ujRs3KiUlpU0+CwCAzsjnkrJlZWWtfvOYmBi5XC5Pu7q6WtHR0ZK+Hb3HxcUpPj5eQUFBGj58\nuL744otWfxYAAGhGsMfHx2vNmjUqKytTeXm551dzjBgxQoWFhZKk0tJSORwOhYWFSZKCg4MVHx+v\n/fv3e4737du3td8DAACoGQvUvPPOOxf0WZYlp9Pp882Tk5OVlJSkjIwMBQcHa/bs2SooKFB4eLjS\n0tKUm5urGTNmyO12q3///hozZkzrvgUAAJAkWe6zH3y3YxUVFUpNTZXT6VR8PDNNAaCl2Jimc/A5\nYq+urtZzzz2n//znP7IsS0OGDNGvfvUrRUVF+aM+AEAbYGOazsPnM/bZs2crKSlJzz77rJ5++mld\neeWVys3N9UdtAIA2crGNaWAWnyP2EydOaMKECZ52//799cEHH9haFACgbbExTefhc8R+4sQJVVdX\ne9qVlZWqr6+3tSgAQNvq7WUDGjamMY/PEfvUqVN15513Kjo6Wm63W1999ZUWLFjgj9oAAG2EjWk6\nD5/BftNNN+n999/X3r17JUl9+/ZVaGio3XUBANoQG9N0Hs2aFf/yyy/riy++kGVZGjBggCZPnqwe\nPXr4oz4AQBthY5rOwecz9pycHF1yySXKysrSxIkTFRQUpIceesgftQEAgBbyOWJ3u9168MEHPe1R\no0Zp0qRJthYFAABax+eI/eqrr9Znn33mae/atUsDBgywtSgAANA6PkfsGzdu1IoVKxQZGanGxkYd\nPXpUDodDhYWFsixLGzZs8EOZAACgOXwG+yuvvOKHMgAAQFvwGeyxsbFat26ddu7cKUkaMmSI0tPT\nbS8MAAC0nM9gnz9/vg4fPqwbbrhBbrdb7777rj755BP95je/8Ud9AGA8dl1DW/IZ7Lt379aKFSs8\n7YkTJ2r8+PG2FgUAnQW7rqGt+ZwVf+rUKTU2Nnrap0+f1unTp20tCgA6C3ZdQ1vzOWJPSUnRXXfd\npeuvv16StHXrVo0dO9b2wgCgM2DXNbS1Zm0Cc+ONN6qkpESWZemJJ57Q4MGD/VEbABivtyNce7+s\nuaCfXdfQWj6DfcGCBZo1a5aGDBnij3oAoFNh1zW0NZ/BHhwcrKKiIg0dOlRdunTx9AcF+Xw8DwDw\ngV3X0NYst9vtvtgLrrvuOn3zzTc6+2WWZZ2zzKw/VFRUKDU1VU6nU/Hx/A8PAEBTfI7Yt2+/8BYR\nAABon9iPHQAAg/gM9pycHF1//fXKysqS2+3W9u3b9dBDD2n58uX+qA8AjMDqcvAX9mMHAJuxuhz8\nif3YAcBmrC4Hf2rWfuwrV65UREQE+7EDQCuwuhz8if3YAcBmrC4Hf/IZ7HFxcf6oAwCMxepy8Cef\nwQ4A+G5YXQ7+RLADgB+MSo4nyOEXFw32Tz75RBs3blR1dbUsy1JsbKxGjx6tpKQkf9UHAABawOuP\nuz3//POaP3++unbtqiFDhujaa6+VJOXm5jKhDgCAdsrriH3Tpk1atWrVOTu6SdK9996rrKwsTZ48\n2e7aAABAC110gZqmtma1LEuNjY22FQQAAFrP64h91KhRGjdunMaMGaPo6GhJ324I8/777+v222/3\nW4EAAKD5Lrof+44dO7Rp0yZVV1dLkq644gqlpKQoMTHRbwWewX7sAAD4dtFZ8TExMXI4HLIsS0FB\nQYqNjfWM3gEAQPvj9Rn76tWrNXHiRO3YsUOnTp1SXV2dtm3bpoyMDK1fv96fNQIAgGbyOmJfvXq1\n3nrrLV122WXn9B85ckS/+MUvlJ6ebntxAACgZbyO2IOCgi4IdUmKiIiQZVm2FgUAAFrH64g9MTFR\n9913n3784x/r8ssvl/TtrPh33nlHQ4cO9VuBAACg+bwG++OPP653331X//znPz2z4mNjY5WRkaG0\ntLRmf8DChQtVUlIiy7KUm5urQYMGXfCaZ555Rp988oleffXVVnwFAABwhtdgtyxLY8eO1dixYy84\nVlRUpOHDh/t88+LiYu3bt0/5+fkqKyvTrFmzlJ+ff85rysrK9NFHH12wwh0AAGi5i648583vf//7\nZr2uqKjIM7pPSEhQTU2Namtrz3lNXl6ecnJyWlMGAAA4j9cR+5IlS5rsd7vdqqioaNabu1wuXXPN\nNZ52ZGSkXC6XunXrJkkqKCjQDTfcoJ49e7akZgAA4IXXEfvbb7+tI0eOKDg4+JxfISEhrZ4Vf/Yi\nd0ePHtXatWs1ZcoUud1uXWQBPAAA0ExeR+yLFy/WkiVLNGfOnAuCfOvWrc1685iYGLlcLk+7urra\ns3Ldli1bdOTIEU2YMEEnT55UeXm58vLyNGPGjNZ8DwAAoIuM2AcPHqy5c+eqvr7+gmNTpkxp1puP\nGDFChYWFkqTS0lI5HA6FhYVJkm6++WatX79e+fn5WrZsmRITEwl1AAC+I68jdrfbrV69eqmxsfGC\nbVpvuummZr15cnKykpKSlJGRoeDgYM2ePVsFBQUKDw9v0Y/MAQCA5vG6u1tWVpaWL1+ugQMHyrKs\nc56BW5alzz77zG9FSuzuBgBAc3gdsS9fvlyStGvXLr8VAwAAvpuLbtsqfTvh7eWXX9YXX3why7I0\nYMAATZ48WT169PBHfQAAoAV8LlCTk5OjSy65RFlZWZo4caKCgoL00EMP+aM2AADQQj5H7G63Ww8+\n+KCnPWrUKE2aNMnWogAAQOv4HLFfffXV50yU27VrlwYMGGBrUQAAoHV8jtg3btyoFStWKDIyUo2N\njTp69KgcDocKCwtlWZY2bNjghzIBAEBz+Az2V155xQ9lAACAtuA12IuLiy/6B6+//vo2LwYAAHw3\nXoM9MzNTV155pQYPHtzkpi8EOwAA7Y/XYF+xYoXWrl2r7du366abbtJtt92mpKQkf9YGAABayGuw\nDxs2TMOGDVNdXZ0KCwu1aNEiuVwupaen69Zbb1VcXJw/6wQAAM3gda34850+fVpr1qzRs88+K6n5\nW7e2FdaKBwDAN5+z4svKyrRmzRq99957SkxM1BNPPKHRo0f7ozYAANBCXoP99ddf19q1a2VZlm67\n7TYVFBQoIiLCn7UBAIAW8norfuDAgerTp49iYmK+feF5M+PP7P7mL9yKBwDAN68jdqfT6c86AABA\nG/Aa7Mx6BwCg4/E5eQ4wzaaPK7TauVv7q46ptyNc41L7aVQyj3cAmIFgR6ey6eMKLVqx3dPe+2WN\np024AzCBz21bAZOsdu5uUT8AdDQEOzqV/VXHmuwv99IPAB0NwY5OpbcjvMn+Xl76AaCjIdjRqYxL\n7deifgDoaJg8h07lzAS51c7dKq86pl7MigdgGIIdnc6o5HiCHICxuBUPAIBBCHYAAAxCsAMAYBCC\nHQAAgzB5DrAJa9IDCASCHbABa9IDCBRuxQM2YE16AIFCsAM2YE16AIFCsAM2YE16AIFCsAM2YE16\nAIHC5DnABqxJDyBQCHbAJqxJDyAQuBUPAIBBCHYAAAxCsAMAYBCCHQAAgxDsAAAYxPZZ8QsXLlRJ\nSYksy1Jubq4GDRrkObZlyxYtXrxYwcHB6tu3rxYsWGB3OQAAGM3WEXtxcbH27dun/Px8zZ8//4Lg\nnjNnjpYuXarXXntNx48f16ZNm+wsBwAA49ka7EVFRUpLS5MkJSQkqKamRrW1tZ7ja9euVUxMjCQp\nKipKX3/9tZ3lAABgPFuD3eVyKSoqytOOjIyUy+XytLt16yZJqq6u1ubNm5WSkmJnOQAAGM+vk+fc\nbvcFfYcPH9b999+vuXPnqnv37v4sBwAA49ga7DExMeeM0KurqxUdHe1pHz9+XPfee69ycnI0fPhw\nO0sBAKBTsDXYR4wYocLCQklSaWmpHA6HwsLCPMfz8vI0ZcoUjRgxws4yAADoNGz9cbfk5GQlJSUp\nIyNDwcHBmj17tgoKChQeHq6RI0fq7bff1v79+/XGG2/IsizdeuutGjdunJ0lAQBgNMvd1IPvdqii\nokKpqalyOp2Kj2fHLAAAmsLKcwAAGIRgBwDAIAQ7AAAGIdgBADAIwQ4AgEEIdgAADEKwAwBgEIId\nAACDEOwAABiEYAcAwCAEOwAABiHYAQAwCMEOAIBBCHYAAAxCsAMAYBCCHQAAgxDsAAAYhGAHAMAg\nBDsAAAYh2AEAMAjBDgCAQQh2AAAMQrADAGAQgh0AAIMQ7AAAGIRgBwDAIAQ7AAAGIdgBADAIwQ4A\ngEEIdgAADEKwAwBgEIIdAACDEOwAABiEYAcAwCAEOwAABiHYAQAwCMEOAIBBCHYAAAxCsAMAYBCC\nHQAAgxDsAAAYhGAHAMAgBDsAAAYJsfsDFi5cqJKSElmWpdzcXA0aNMhzbPPmzVq8eLGCg4M1atQo\nTZ061e5yAAAwmq3BXlxcrH379ik/P19lZWWaNWuW8vPzPccXLFigl156STExMZo4caJuvvlmJSQk\n2FkSWmHTxxVa7dyt/VXH1NsRrnGp/TQqOT7QZQEAmmDrrfiioiKlpaVJkhISElRTU6Pa2lpJUnl5\nuSIiIuRwOGRZllJSUrRlyxY7y0ErbPq4QotWbNfeL2vU2OjW3i9rtGjFdm36uCLQpQEAmmBrsLtc\nLkVFRXnakZGRcrlcTR6LiopSdXW1neWgFVY7d7eoHwAQWH6dPOd2u1t1DIGzv+pYk/3lXvoBAIFl\na7DHxMR4RuiSVF1drejoaM+xQ4cOeY5VVVUpJibGznLQCr0d4U329/LSDwAILFuDfcSIESosLJQk\nlZaWyuFwKCwsTJIUFxen2tpaHTx4UA0NDdqwYYNGjhxpZzlohXGp/VrUDwAILFtnxScnJyspKUkZ\nGRkKDg7W7NmzVVBQoPDwcKWlpWnOnDnKycmRJKWnp6tPnz52loNWODP7fbVzt8qrjqkXs+IBoF2z\n3B3k4XZFRYVSU1PldDoVH0+oAADQFFaeAwDAIAQ7AAAGIdgBADAIwQ4AgEEIdgAADEKwAwBgEIId\nAACDEOwAABiEYAcAwCAEOwAABiHYAQAwCMEOAIBBCHYAAAxCsAMAYBCCHQAAgxDsAAAYhGAHAMAg\nBDsAAAYh2AEAMAjBDgCAQUICXUBznT59WpJUWVkZ4EoAAPCv2NhYhYQ0L7I7TLAfOnRIkjRhwoQA\nVwIAgH85nU7Fx8c367WW2+1221xPm6irq9POnTsVHR2t4ODgQJcDAIDftGTE3mGCHQAA+MbkOQAA\nDEKwAwBgEIIdAACDEOwAABik3f64W11dnWbMmKHDhw+rvr5e999/vwYOHKhHH31Ubrdb0dHReuqp\np9SlS5dAl4omNHX+CgsLtXPnTkVGRkqSsrOzlZKSEuBK4c3JkyeVnp6uBx54QN///ve59jqYs8/f\n1q1bufY6gG3btunBBx9Uv3795Ha7NWDAAN1zzz0tvvbabbB/8MEHGjRokLKzs3Xw4EFNmTJFQ4cO\n1cSJE3XzzTdr8eLFevPNN5WRkRHoUtEEb+fvkUce4S+UDuKFF15QRESEJGnJkiXKzMzUj370I669\nDuLs8yeJa6+D+N73vqclS5Z42jNnzmzxtddub8WPHTtW2dnZkqSDBw/qiiuuUHFxscaMGSNJGj16\ntDZv3hzIEnERTZ0/SeKnKzuGPXv2aM+ePUpJSZHb7VZxcbFGjx4tiWuvIzj//Elcex3F+edp27Zt\nLb722m2wn5GRkaHHHntMM2fO1IkTJzy3IHr06OFZjQ7t15nzl5ubK0lauXKlJk2apIcfflhff/11\ngKuDN08++aRmzJjhaXPtdSxnnz/LsiRx7XUUZWVlmjp1qiZMmKDNmzerrq6uxddeu70Vf0Z+fr52\n7dqlRx555Jx/yfCvz47h7POXm5uriIgIDRw4UH/84x+1dOlSPf7444EuEed56623lJycrLi4uCaP\nc+21b+efP7fbrdtvv51rrwPo06ePpk2bpltuuUXl5eXKyspSQ0OD53hzr712G+ylpaXq0aOHYmNj\nNXDgQDU2Nqpbt26qr69X165dVVVVpZiYmECXCS/OP3+nT59W//79FRUVJUlKTU3V3LlzA1skmrRx\n40ZVVFToH//4h6qqqtSlSxeFhYVx7XUQZ5+/yspKhYaG6v/+7/80cOBASVx77ZnD4dAtt9wiSerV\nq5cuv/xy7dy5s8XXXru9FV9cXKyXXnpJkuRyufTNN99o+PDheu+99yRJhYWF+sEPfhDIEnERTZ2/\nOXPmqLy8XJK0detW9e/fP5AlwovFixdr9erVev3113XXXXfpgQce4NrrQM4+f+PGjdPUqVO1atUq\nrr0OYN26dZ6/Nw8dOqTDhw/rzjvvbPG1127Xij958qRyc3NVWVmpkydPavr06UpKStJjjz2m+vp6\n9ezZUwsXLmRDmHbq/PM3bdo0hYWF6amnntKll16qbt266be//a1nBI/2admyZYqPj9fIkSO59jqg\nM+evZ8+eXHsdQG1trR5++GEdO3ZMDQ0NmjZtmgYOHKhf//rXLbr22m2wAwCAlmu3t+IBAEDLEewA\nABiEYAcAwCAEOwAABiHYAQAwCMEOAIBBCHbAYAcOHGjRjl4bN25UTU2NJOnhhx9WdXW1JOntt99u\n0ec+99xzWrZsWYv+DIC2QbADhjuzCUhz/PnPf/ZsEPLMM88oJiZGVVVVys/Pt6s8AG2MYAc6od/9\n7ncaP368xo8fr5ycHJ0+fVqrVq3SRx99pEcffVRlZWUaM2aMysvL9cgjj2j37t2aMWOGtm3bpvHj\nx3veZ+bMmVqzZo2kb5cyvf3223Xfffdp3759ntds2bJFmZmZyszM1M9+9jMdOHDA798X6EwIdqCT\nOX36tC699FK99tpreu2113T06FF9+OGHuvvuu3X55Zfr6aefVkJCgmekP336dPXv3195eXmSmr4D\nsHfvXq1fv15vvvmmnn/+eU+w19XVae7cuXr++ef16quvasKECXryySf992WBTqjd7u4GwB7BwcEK\nCgrS3XffrZCQEP3vf//TkSNHLnhdS1ab/vzzz5WUlKSQkG//Shk2bJin/9ChQ5o2bZrcbrfcbneL\nHg0AaDmCHehk/v3vf2vt2rVau3atQkND9ctf/rLZf/b8UK6vr5f07T8CgoL+/w3AxsZGSVLXrl3V\ns2dPLV++vA0qB9Ac3IoHDHf+yPvw4cOKi4tTaGioDhw4oJKSEk9ABwUF6dSpU+e8/uy+yy67TFVV\nVZKkEydOaMeOHZKkhIQEffrpp2poaNCpU6e0bds2SVLfvn115MgR7d69W9K32/m+8cYb9n1ZAIzY\nAdMdOXJEWVlZntvgSUlJOn78uCZMmKCrrrpK06dP1wsvvKAbbrhBI0eO1P3336+8vDzP6Pyqq66S\ny+VSdna2XnzxRfXv31933nmnevfuraFDh3pek5qaqp/+9Kfq2bOnEhMTJUmhoaFatGiRZs2apdDQ\nUEnSvHnzAvMfAugk2LYVAACDcCseAACDEOwAABiEYAcAwCAEOwAABiHYAQAwCMEOAIBBCHYAAAxC\nsAMAYJD/Byq/ngA5VRVBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8812f74c88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.sort_values('Latitude').plot('Latitude', 'p, Mpi100', ls='', marker='o')\n",
    "plt.ylabel('Mpi100 proportion')\n",
    "plt.legend().set_visible(False)\n",
    "plt.xlim(30, 50)\n",
    "plt.ylim(-0.1, 1.1)\n",
    "sns.despine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Allele</th>\n",
       "      <th>Latitude</th>\n",
       "      <th>Location</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>48.1</td>\n",
       "      <td>Port Townsend, WA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>48.1</td>\n",
       "      <td>Port Townsend, WA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>48.1</td>\n",
       "      <td>Port Townsend, WA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>48.1</td>\n",
       "      <td>Port Townsend, WA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>48.1</td>\n",
       "      <td>Port Townsend, WA</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Allele  Latitude           Location\n",
       "0       0      48.1  Port Townsend, WA\n",
       "1       0      48.1  Port Townsend, WA\n",
       "2       0      48.1  Port Townsend, WA\n",
       "3       0      48.1  Port Townsend, WA\n",
       "4       0      48.1  Port Townsend, WA"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rows = []\n",
    "for i, row in df.iterrows():\n",
    "    for _ in range(row['Mpi90']):\n",
    "        rows.append({'Location':row['Location'], 'Latitude': row['Latitude'], 'Allele': 0})\n",
    "    for _ in range(row['Mpi100']):\n",
    "        rows.append({'Location':row['Location'], 'Latitude': row['Latitude'], 'Allele': 1})\n",
    "raw_df = pd.DataFrame(rows)\n",
    "raw_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-7.59881547] [[ 0.17754715]]\n"
     ]
    }
   ],
   "source": [
    "clf = sklearn.linear_model.LogisticRegression(C=1e12, random_state=0)\n",
    "clf.fit(raw_df['Latitude'].reshape(-1, 1), raw_df['Allele'])\n",
    "print(clf.intercept_, clf.coef_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is very close to McDonald's intercept of -7.6469 and slope of 0.1786."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7f8837430be0>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGBCAYAAACTuDAhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0VeW9//HPyclAEjJC5pOEEIZASGQOISiVRHHqqvRe\nvRRrvWr7W3WqFq0VWFV7K0rba9Wrda32tt7eutRcW8FVpSUKKmoGCCBgwmjmk/mEkJB5Or8/kFMw\nAwSTc7KT92stFuz9POfsbzjs5MPez34ek91utwsAAMCg3FxdAAAAwNdBmAEAAIZGmAEAAIZGmAEA\nAIZGmAEAAIZGmAEAAIbmsjBz4sQJXXPNNXrttdf6teXl5enf/u3ftG7dOm3atMkF1QEAAKNwSZhp\nb2/XU089pdTU1AHbn3jiCb344ot6/fXX1dLSoo8//tjJFQIAAKNwSZjx8vLSH/7wB4WGhg7YvnXr\nVkdbcHCwTp8+7czyAACAgbgkzLi5ucnT03PQdl9fX0lSXV2dcnJytHLlSmeVBgAADGbMDgBuaGjQ\nPffcoyeffFIBAQGD9uvp6ZHValVPT48TqwMAAGPFmAwzLS0t+sEPfqD169cPOq7mnJqaGqWnp6um\npsZJ1QEAgLFkTIaZLVu26M4771RaWpqrSwEAAGOcuysOWlhYqC1btqiqqkru7u7KysrSqlWrZLFY\ntGLFCv3tb39TeXm53nzzTZlMJn3zm9/ULbfc4opSAQDAGOeSMJOYmKhXX3110PbDhw87sRoAAGBk\nY/I2EwAAwKUizAAAAEMjzAAAAEMjzAAAAEMjzAAAAEMjzAAAMEFUVlYqISGh31PD//Iv/6INGzYM\n+rpPPvlEmZmZkqSmpiZ9//vf14MPPuho7+np0SOPPKJ169bp9ttvl9VqlSQdO3ZMa9eu1bp16/Tz\nn/98FL6iswgzAABMIDExMXr33Xcd2+Xl5Tpz5syQr7nyyiu1du1aSdKTTz6pxYsXX9D+7rvvKiAg\nQK+//rp++MMf6tlnn5UkPf300/rZz36m119/Xc3Nzfrkk09G+Ks5yyXzzAAAMFFll+dr25EsWZur\nZfGP0Jq5q5UWs+Sy32/btm365JNP1NLSotraWt1xxx369re/PWj/5ORk5eTkyG63y2Qyafv27Vqx\nYoXa29slSatWrdKaNWuUl5cnT09Pvfjii3r//fd14sQJ/fSnP9XmzZtVUFCgo0ePOt4zNzdXN998\nsyRp+fLl2rRpk7q7u2W1WpWYmOh435ycHF155ZWX/bUOhiszAAA4SXZ5vl7IfUXlTZXqs/epvKlS\nL+S+ouzy/K/1vl988YV+97vf6U9/+pOef/75Ift6eHgoOTlZeXl5kqRdu3Zp5cqVF/SZMWOGXnvt\nNSUkJGjr1q2SJJPJJEny8fHp9542m03BwcGOfiaTSTabTYGBgY4+wcHBqq+vv/wvcgiEGQAAnGTb\nkawB9789yP5LtXTpUplMJgUFBSkgIECnTp0asv91112nd999VydPnlR4eLi8vb0vaD+3yPP8+fNV\nWlo67HrOXfWx2+3Dfu3lIMwAAOAk1ubqYe2/VH19fY4/nwsSQ0lNTdWePXu0fft2rV69etD3u5T3\nkqTQ0FDZbDZJZwcD2+12hYSE6PTp044+tbW1Cg0NvaSvZ7gIMwAAOInFP2JY+y/VwYMHZbfbderU\nKbW1tSkoKGjI/h4eHlqyZIneeustXX311f3a9+/f73jf+Pj4fu12u/2Cqy5paWnasWOHJOmDDz5Q\nSkqKzGazpk+frgMHDkiS3nvvvVEZLyMxABgAAKdZM3e1Xsh9pd/+m+f2vzoyHFFRUfrRj36k8vJy\n/fjHP5Yk/f73v1dKSoquuOKKAV9z3XXXqbGxUZMnT+7XVlBQoNdee01ubm760Y9+pKyss7fB+vr6\ndMcddzgGG3/ve9/TfffdpxtuuEHZ2dlat26dvLy8tGXLFknSxo0b9fjjj8tut+uKK65w3L4aaSa7\ns25ojRKr1ar09HTt2rVLFovF1eUAADCk7PJ8vX3e00w3j8DTTCdPntSjjz56wf7du3fLx8dHS5YM\n771XrVql7du39xtHM5ZxZQYAACdKi1nytcLLpfL09Bz0qsxQLmWMzFjDlRkAAGBoDAAGAACGRpgB\nAACGRpgBAACGRpgBAACGRpgBAACGRpgBAACGRpgBAACGRpgBAACGRpgBAACGRpgBAACGRpgBAACG\nRpgBAACGRpgBAACGRpgBAACGRpgBAACGRpgBAACGRpgBAACGRpgBAACGRpgBAACGRpgBAACGRpgB\nAACGRpgBAACGRpgBAACG5rIwc+LECV1zzTV67bXX+rXl5OTolltu0dq1a/Xyyy+7oDoAAGAULgkz\n7e3teuqpp5Samjpg++bNm/XSSy/pjTfeUHZ2toqKipxcIQAAMAqXhBkvLy/94Q9/UGhoaL+2iooK\nBQYGKiwsTCaTSStXrlReXp4LqgQAAEbg7oqDurm5ydPTc8A2m82m4OBgx3ZwcLAqKiqcVRowInbu\n3Cmr1SqLxaKMjAynHTcvL0/V1dWKiIjQsmXLnHZcGJ/JZHL82W63D7pvqP7O4KrjGslE/Dsa8wOA\nJ8oHgfFj586dam1tVVBQkFpbW7Vz506nHDcvL0+tra3y9/dXa2srVzRxyc7/4Xdue6B9Q/V3Blcd\n10gm6t/RmAszoaGhqq+vd2zX1tYOeDsKGKusVuuQ26Olurp6yG0AGK/GXJiJiopSa2urqqqq1NPT\no48++kgrVqxwdVnAJbNYLENuj5aIiIghtwFgvHLJmJnCwkJt2bJFVVVVcnd3V1ZWllatWuUYX/DE\nE09o/fr1kqSbbrpJsbGxrigTuCwZGRkuGTOzbNkyxszgstjt9mGNmRms/2hz1XGNZKL+HZnsBv9K\nrVar0tPTtWvXLqf9DxgAAIwdY+42EwAAwHAQZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKER\nZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAA\ngKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKER\nZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAAgKERZgAA\ngKERZgAAgKERZgAAgKERZgAAgKERZgAAgKG5u+KgzzzzjA4dOiSTyaSNGzcqKSnJ0fbaa6/pnXfe\nkdls1rx587RhwwZXlAgAAAzC6WEmPz9fZWVlyszMVFFRkTZt2qTMzExJUktLi/74xz9q165dMplM\nuvvuu3X48GElJyc7u0wAAGAQTr/NlJubq4yMDElSfHy8mpub1draKkny9PSUp6enWlpa1NPTo46O\nDgUEBDi7RAAAYCBODzM2m03BwcGO7aCgINlsNklnw8x9992njIwMpaenKzk5WbGxsc4uEQAAGIjL\nBwDb7XbHn1taWvS73/1O7733nnbt2qVDhw7p+PHjLqwOAACMdU4PM6GhoY4rMZJUV1enkJAQSVJx\ncbGio6MVEBAgd3d3LV68WIWFhc4uEQAAGIjTw0xaWpqysrIkSYWFhQoLC5OPj48kKSoqSsXFxerq\n6pIkFRQUcJsJAAAMyelPMy1YsECJiYlau3atzGazHn/8cW3btk1+fn7KyMjQ3Xffrdtvv13u7u5a\nsGCBFi1a5OwSAQCAE/X22XWs9JT2H6tVQVGDfvXAlcN6vcl+/qAVA7JarUpPT9euXbtksVhcXQ4A\nALgELe3d+uxYnfYerdH+o7U609btaHvn2W8N671cMmkeAACYeCrrW5R/pEb5R2pVWNyg3r6RuZ5C\nmAEAAKOip7dPR0oalH+kVvlHalRZ3zpo32B/Ly2eE67Fc0KHfRzCDAAAGDHNrV3af6xW+UdqdeBY\nrVo7egbtOyM6UEvnhGlJYrjiowJkMpku65iEGQAA8LXUNLQqr6BaeQU1OlrSoMHuHnl5mjV/ZoiW\nJoZr8ZwwBftPGpHjE2YAAMCw2O12lVY3K/fzauUVVKukqnnQviFB3loyJ0xLE8OVFD9Vnh7mEa+H\nMAMAwDiRXZ6vbUeyZG2ulsU/QmvmrlZazJIRee9zj0+fCzC1p9oG7GcySbNjgrRkbriWJoYrNtzv\nsm8fXSrCDAAA40B2eb5eyH3FsV3eVOnYvtxA09Xdq0Mn65X7ebX2HqlRU0vXgP3czW6aPytEqUkR\nWjo3XIF+Xpd1vMtFmAEAYBzYdiRrwP1vH8kaVphpbe/WvqO1yi2o1oFjtWrv7B2wn88kdy2eE6bU\npAgtnB0qn0kel1X3SCDMAAAwDlibq4e1/3wt7d3aU1CtTw9V6eCJOvX0DjyCN8jPSynzIpQ6L0JJ\nM6bKw93l61VLIswAADAuWPwjVN5UOeD+gbS0dSmvoEbZh4cOMBFTfZU6L0KpSRGaFRMkN7fRHf9y\nOQgzAACMA2vmrr5gzMw5N89d7fjz2QBz9grMoZP1gwaY6VEBWp4UoWVJEYoJG/0BvF8XYQYAgHHg\n3LiYt897munmuauVPPUKvb+nTJ8ertKhE/WDLiEwwxKgtCuilJYcqYipvs4s/WsjzAAAME6kxSxR\nWswSNbeevQLz3o4q/frkjsEDTHSgViRHKu2KSIVPMVaAOR9hBgCAcaCto1t5BTX6+DOrDg5xBWZW\nTKDSkqOUdkWkwoJ9nFzl6CDMAABgUF3dvdp3tFYff1ap/CM16urpG7Df7JggpV0RqbTkSIWOkwBz\nPsIMAAAG0tvbp0Mnbdr9mVV5BdVqG2Qhx9kxQVoxP1LLkyMVGjT+Asz5CDMAAIxxfX12HS09pY8/\nsyr7cNWgM/FOi/DXVQuidOX8KEOPgRkuwgwAAGOQ3W5XcWWTPv6sUp8cqlR9Y/uA/cKn+GjlAouu\nWhClmHB/J1c5NhBmAAAYQ2pPtemjAxX6aL9V1rqWAfsE+3tpxfworVxg0czowDE/D8xoI8wAAOBi\nLe3dyj5UqQ/3W1VY3DBgHz8fDy1PjtRVC6KUOH2qzGNwJl5XIcwAAOAC3T19OnCsVh/ut2rvkRp1\nD/Ak0iRPs1ISI7RyYZTmzwodM2shjTWEGQAAnMRut+tEeaM+3G/Vx59V6kxb/4G8biZp/uxQXb0o\nWssSwzXJix/VF8PfEAAAo6ymoVUf7rfqo/0VqrK1Dtgn3hKgqxdF66r5UQryn+TkCo2NMAMAwCho\nae/WJwcr9eG+Ch0tPTVgn6mB3vrGQouuXmSZsE8ijQTCDAAAI6S3z65DJ+u1a2+5cguqBxwH4+3l\nrrTkSF292KJ506fKjYG8XxthBgCAr6mqvkU788v14b4K2Zo6+rW7uZm0cHaoVi2K1tJ54fLyMLug\nyvGLMAMAwGVo6+jWp4eqtHNv+aC3keItAVq1OFpXzbco0M/LyRVOHIQZAAAuUV+fXQXFNu3cW66c\nz6vV2dXbr4+/r6e+sciijCUxiosMcEGVEw9hBgCAi6hpaNUH+yq0a1+F6k619Wt3czNpyZwwpS+J\n0eI5YcwH42SEGQAABtDR2aOcz6u0c2+FPi+yDdgnNtxPGUtjtHKhRUF+PE7tKoQZAAC+ZLfbdbLi\ntN7bU6aPP7OqvbP/baTJ3h5aufDsbaR4S8CEXxdpLCDMAAAmvDNtXfpov1Xv7SlTaXVzv3Y3k7Qw\nIUzpS6KVkhguD3eeRhpLCDMAgAnJbreroKhBWXllyvm8asA5YaJCJitjaYyuXmTRlABvF1SJS0GY\nAYAJLrs8X9uOZMnaXC2Lf4TWzF2ttJglri5r1Jxq7tCu/HK9v7dc1QMsLeDpYdaV8yN1bUqs5kwL\n5jaSARBmAGACyy7P1wu5rzi2y5sqHdvjKdD09vbpwPE6ZeWVKf9orfr67P36zLAE6NqUWF21wCJf\nbw8XVInLRZgBgAls25GsAfe/fSRrXISZmoZW7dxbrp355WoYYGZe30nu+saiaF2zNEbxlkAXVIiR\nQJgBgAnM2lw9rP1G0N3Tq7yCGr2XV6aDJ+sH7JM4fYquTYnV8uQITfLkR6HR8QkCwARm8Y9QeVPl\ngPuNprK+RTtyS7Urv0Jn2rr6tQdO9lL6kmhdkxKrqJDJzi8Qo8YlYeaZZ57RoUOHZDKZtHHjRiUl\nJTnaampqtH79evX09Gju3Ll68sknXVEiAEwIa+auvmDMzDk3z13tgmqGr7unT3kF1dqRW6rDX/Sf\n2M5kkhbODtXqZbFaMjdc7mZm5h2PnB5m8vPzVVZWpszMTBUVFWnTpk3KzMx0tG/ZskV333230tPT\n9Ytf/EI1NTUKDw93dpkAMCGcGxfz9nlPM91sgKeZahpalZVXpp17y3W6pbNfe0iQt65ZGquMJTEK\nCeKR6vHO6WEmNzdXGRkZkqT4+Hg1NzertbVVvr6+stvt2r9/v5577jlJ0s9+9jNnlwcAE05azJIx\nH16ks08k7T1Sqx25pfrsRJ3sX3kgyc0kLZkbruuXT9P8WaEyu/FI9UTh9DBjs9k0b948x3ZQUJBs\nNpt8fX116tQp+fj4aPPmzTpy5IgWL16s9evXO7tEAMAYUt/Yrvf2lOm9PWU61dz/iaQpAZO0OiVW\n16TEamogV2EmIpcPALafF63tdrvq6ur07//+74qMjNT/+3//T7t379bKlStdWCEAwNl6++z67Hid\nduSWKv9Ijb46Lcy5sTDXp07T4jlhMjMWZkJzepgJDQ2VzfbPQVp1dXUKCQmRdPYqTVRUlCwWiyQp\nNTVVX3zxBWEGACaIU80den9vmd7LK1NdY3u/9kA/L12zNEarl01TWLCPCyrEWOT0MJOWlqaXXnpJ\nt956qwoLCxUWFiYfn7P/IM1msywWi8rLyxUTE6PCwkLddNNNzi4RAOBEfX12Hf6iXv/ILdWeghr1\nDjA77xUzp+r61DgtTQyXhztXYXAhp4eZBQsWKDExUWvXrpXZbNbjjz+ubdu2yc/PTxkZGdq4caMe\ne+wx2e12zZo1S6tWrXJ2iQAAJ2hq6dSu/HLtyCsbcI0kPx9PZSyN0XXLYhXJvDAYgslu/+p4cGOx\nWq1KT0/Xrl27HLenAABjk91u1/GyRm3PLtGnh6rU09t/perE6VN0Xeo0LU+KkKeH2bF/oi2IiUvn\n8gHAAIDxr7O7Vx8fsGp7TomKrE392n0nuWvVkhitXhar2HD/fu0TZUFMXB7CDABg1FTbWvX3nBLt\n3Fuulvbufu2zYgJ1fWqcVsyPHHKNpPG+ICa+HsIMAGBE9fXZdeB4nbZnl2j/sdp+k9t5urvpqgUW\n3ZgWpxnRl7ZS9XhcEBMjhzADABgRZ9q69P6ecv0jt0Q1DW392sOCfXTD8jhlLI2Rv6/nsN57PC2I\niZFHmAEAfC1fWE/r79kl2n3Aqq6eCwf0npvc7sa0OC1MCLvsJQaMviAmRhdhBgAwbN09vco+VKV3\ns0t0vKyxX/tkbw9lLI3RDcvjFDHV92sfz6gLYsI5CDMAgEtW19imHbmlem9PmZpauvq1T48K0I1p\ncbpqQdSQA3ovh1EWxITzEWYAAEOy2+06fNKmd7OLtbew/zpJ7maT0pKjdNOKOM2ODZLJxGrVcC7C\nDABgQG0d3dqVX6G/55TIWtfSr31qwCRdt3yark2JVZDfJBdUCJx10TDT1dWlv/zlL6qurtYjjzyi\nQ4cOKSEhQV5eXs6oDwDgZGU1zdqeXaIP91Woo6u3X/sVM6fqxrQ4LZ0bzmrVGBMuGmaefPJJ+fn5\n6cCBA5KkwsJC/elPf9Jzzz036sUBAJyjp7dPeQXV2p5dooKihn7t3l7uSl8crRvS4hQd5ueCCoHB\nXTTMFBcXKzMzU7fffrskad26ddq+ffuoFwYAGH2nmjuUlVuqHXllOtXc0a89OsxPN62I0zcWWuQz\nycP5BQKX4KJhxt39bJdzA7ra2trU0dH/HzwAwBjsdruOlJzS9uwS5RyuUu9XRvS6uZmUOi9CN6bF\naV78lH4DelnwEWPNRcPMddddpzvuuENWq1VPPfWUPv74Y61bt84ZtQEARlBHZ48+OmDV9uwSlVY3\n92sP9PPSdcum6brUWE0J8B7wPVjwEWPRRcPMd7/7XSUnJ2vv3r3y9PTUb37zG82bN88ZtQEARkBl\nfYv+nl2iXfnlau3o6dc+Ny5YN6bFKTUpUh7uQw/oZcFHjEWDhpnc3NwLthMTEyVJZ86cUW5urlJT\nU0e3MgDAZevts2vfkRq9m12igyfq+7V7eZr1jYVnF3uMiwy45PdlwUeMRYOGmZdffnnQF5lMJsIM\nAIxBTS2dem9PmXbklqqusb1fe+RUX92QFqf0JTGa7D38Ab0s+IixaNAw8+qrr16wbbfbmdURAMao\nE+WN2p5dok8OVqp7gMUel8wJ140r4jR/ZojcLnOxR4kFHzE2XXTMzLFjx7Rx40a1tbVpx44d+u1v\nf6sVK1boiiuucEZ9AIBBdHb36pPPKrU9p0RfVJzu1+7n46lrU2J0/fI4hQX7jMgxWfARY9FFw8x/\n/Md/6Omnn9bmzZslSTfccIM2bNigzMzMUS8OANBfTUOrY7HHM23d/dpnRgfqxrQ4XTk/Sp4e5hE/\nPgs+Yqy5pHlmEhISHNtxcXGOuWcAAM7R12fXZyfqtD27RPuO1sr+lcUePdzddOX8KN2YFqdZMUGu\nKRJwkUsKMxUVFY7xMrt375b9q2cRAGBUtLR1aeeXiz1W21r7tYcGeeuG5XHKWBqjgMkjt2YeE+PB\nSC4aZh599FHde++9Kikp0aJFixQVFaVf/vKXzqgNACas4sombc8u0UcHrOrq7r/Y48LZoboxLU6L\n5oTJ/DUG9A6EifFgNBcNMwkJCXrnnXd06tQpeXp6avLkyc6oCwAmnO6ePuUcrtL27BIdLT3Vr913\nkrvSl8bohuVxigoZve/FTIwHoxk0zPzkJz8Z8lHsX/3qV6NSEABMNLbT7dqRW6qsPWU6faazX/u0\nCH/dtCJOKxdYNMlr9McsMjEejGbQs2L58uXOrAMAJhS73a7DX9j095wS5RXUqO8riz2a3UxKS47U\nDWlxmhsX7NR5vpgYD0YzaJhZvHixM+sAgAmhraNbH+w7O6C3oralX/uUgEm6LnWaVqfEKsh/kgsq\nZGI8GM+gYeaOO+6QyWS64Mmlc/8zMJlM2rlz5+hXBwDjRFl1s7Znl+jD/RXq6Oo/oDcpfqpuTItT\nyrxwuZuHXuxxtDExHoxm0DDzwQcfXLDd0dGhrKwsbd26VUVFRaNeGAAYXXdPn/I+r9b2nBIVFjf0\na/f2MuvqRdG6IS1OseH+LqhwcEyMByO56EiygwcPauvWrfrHP/6h3t5e/eIXv9C1117rjNoAwJAa\nmtq1I7dMWXmlahxgQG90mJ9uTIvT1Yss8pk0/MUeAVxo0DDz3//939q2bZva29v1rW99S3/961/1\n4IMP6sYbb3RmfQBgCHa7XZ8X2bQ9e+ABvW5uJqUmRejG5XGaFz+FhXuBETRomHn++ec1Y8YMPf74\n41q2bJkkcfIBwFe0dXTrw30V2p5TqoraM/3ag/29tHrZNK1eFqspAd4uqBAY/wYNMx999JG2bdum\nJ554Qn19fVqzZo26u/svaAYAE1FZzdkBvR/tr1B7Z/8BvfPip+jGtDgtmxfh8gG9wHhnsl/CQkv5\n+fl66623lJWVpZSUFH3nO9/RypUrnVHfRVmtVqWnp2vXrl2yWCyuLgfAONbT26fcz6v195wSFRQN\nPKD3G4uidePyOMVGjK0BvcB4dklh5pyWlha9++672rp1q958883RrOuSEWYAjLaGpnZl5Z0d0Huq\neaABvZN14/I4Xb04mgG9gAsMa17syZMna+3atVq7du1o1QMAY4LdbldBUYO2Z5cot6B64AG98yJ0\nYxoDegFXG/1FPgDAQFrbu/XR/gr9PbdU5TX9B/QG+Z0d0HtdKgN6gbGCMAMAkr6oOK1/5JZq92dW\ndQ4wQ2/i9LMDelOTGNALjDWEGQATVkdnjz45WKl/5JbqZMXpfu2TPM26enG0blgep2kM6AXGLJeE\nmWeeeUaHDh2SyWTSxo0blZSU1K/Ps88+q4MHD+rVV191QYUAxrPymmb9I7dUH+6rUGtHT7/22HA/\nXb+cGXoBo3B6mMnPz1dZWZkyMzNVVFSkTZs2KTMz84I+RUVF2rdvnzw8+CYCYGR09/Qq53C1/pFb\nOuA6SR7ublpxRaSuT41TwrQgBvQCBuL0MJObm6uMjAxJUnx8vJqbm9Xa2ipfX19Hny1btmj9+vV6\n8cUXnV0egHGmpqFVO3JLtTO/XE0tXf3aI6b66vrUaVq1OFoBk72cXyCAr83pYcZms2nevHmO7aCg\nINlsNkeY2bZtm1JSUhQZGens0gCME729fdp7pFY7ckt14Hhdv3Y3N5OWzQvX9anTlDwjRG5uXIUB\njMzlA4DPn7OvqalJW7du1Z/+9CdVV1drGPP5AYAamtr1Xl6ZsvaUqaGpo1/71IBJWp06TdcsjeGx\namAccXqYCQ0Nlc1mc2zX1dUpJCREkpSXl6fGxkbddttt6uzsVEVFhbZs2aLHHnvM2WUCMIjePrsO\nnqhTVl6Z9hT2X63aZJIWJYTp+tRpWpQQKjOPVQPjjtPDTFpaml566SXdeuutKiwsVFhYmHx8fCRJ\nq1ev1urVqyVJlZWV2rBhA0EGwIBsp9u1M79c7+0pU31je7/2wMleuiYlRtemxCp8iu8A7wBgvHB6\nmFmwYIESExO1du1amc1mPf7449q2bZv8/PwcA4MBYCC9vX3ad7RWWXvKtP9orfoGuBOdFD9V16dO\n07KkCHm4cxUGmAiGtdDkWMRCk8D4V3uqTe/vKdP7e8t1qrn/WBg/H0+tWhyt1ctiFR3m54IKAbiS\nywcAA8BAunv6tLewRll5pTp4sl4D/bcrecZUrV4Wq9SkCHm4m51fJIAxgTADYEypqm/Re3vKtCu/\nQqdbOvu1B/p5KX1xtK5dFqvIqZNdUCGAsYYwA8Dlurp7lft5tbLyyvR5ka1fu8kkLZgdqtUpsVqa\nGM5CjwBN/0TVAAAah0lEQVQuQJgB4DLlNc3K2lOmD/dV6Exbd7/2KQGTlLE0RtcsjVVYsI8LKgRg\nBIQZAE7V2t6tTw5Waufech0vb+zX7maSFs8J1+rUWC2azbwwAC6OMANg1PX12VVQbNP7e8uVc7ha\nXd29/fqEBnnr2pRYZTA7L4BhIswAGDV1jW36YF+Fdu4tV+2ptn7t7maTliaGa/WyaZo/kzWSAFwe\nwgyAEdXV3au8gmq9v7dchwZ5pHpahL+uSYnRygUWVqoG8LURZgB8bXa7XUXWJu3ML9dHB6xqbe8/\nmNfX20PfWGhRxtIYxUcFyGTiKgyAkUGYAXDZmlo6tfuAVe/vLVdpdXO/dpNJmj8zRNcsjVXKvHB5\nejCxHYCRR5gBMCw9vX3af7RWH+yv0N7CGvX09r+PFBbso4ylMVq1OFqhQTxSDWB0EWYAXNS520gf\n7K/Q7gNWNbd29evj6WFWWnKErlkaq8TpUxjMC8BpCDMABmU73a6PDlj1wb4KVdSeGbDP7NggZSyJ\n0ZXzo+Tr7eHkCgGAMAPgKzo6e5RbUK0P9lUM+jTS1IBJunpxtK5eFM0q1QBcjjADwDGp3a78CuUc\nrlJHV/9J7SZ5mrU8OVKrFkVr3oypMnMbCcAYQZgBJrCK2jP6cH+FPtxvle10e792k0lKnjFVqxbH\nKDUpQt5efMsAMPbwnQkYYdnl+dp2JEvW5mpZ/CO0Zu5qpcUscXVZDg1N7fr4s0rt/syqImvTgH0s\noZO1anG0vrEwWiFBLC0AYGwjzAAjKLs8Xy/kvuLYLm+qdGy7MtC0tHUp+3C1Pv7Mqs+LbAOOg/Hz\n8dTKBVG6enG0ZkYHMqkdAMMgzAAjaNuRrAH3v30ky+lhprO7V/uO1OqjAxXad7ROPb19/fq4m920\nZG6YVi2O1qKEMHm4s0I1AOMhzAAjyNpcPaz9I623z67DJ+v10QGrcj+vVntnT78+58bBrFxgUWpy\npCbzODUAgyPMACPI4h+h8qbKAfePFrvdrpMVp7X7gFUfH6zU6TOdA/abER2olQssunJ+pKYEMA4G\nwPhBmAFG0Jq5qy8YM3POzXNXj+hx7Ha7Squb9cnBSn16qErVttYB+0VM9dU3Flp01YIoWUKZDwbA\n+ESYAUbQuXExb5/3NNPNI/g0U1nNlwHmYJUq61sG7BPk56UrF0Rp5QILA3kBTAiEGWCEpcUsGdHB\nvhW1Z/TpoSp9crBy0CUFvL3ctTw5Qt9YaFHSjBAmtAMwoRBmgDGoytbiuAJTWt08YJ9JnmYtTQzX\nlfOjtHB2qDw9zE6uEgDGBsIMMEbUNLQ6rsAUVw48mZ2Xp1lL5oRpxfwoLZ4TJi8CDAAQZgBXqqxv\nUc7hKuV8Xq0vKk4P2MfT3U2L5oTpyvlRWjInTJMGWVJgrM88DACjhTADONG5p5ByDlcr5/MqldcM\nPAbG3eymRQmhWjE/Skvnhsln0tBzwYzVmYcBwBkIM8Ao6+uz60R5o3I+r1bu51WqaWgbsJ+72aT5\ns0J15fwopSSGy3cYk9mNpZmHAcDZCDPAKOjt7VNBcYNyP69W7ufVOtXcMWA/T3c3LUwIVWpSpJbO\nDdNkH8/LOp6rZx4GAFcizAAjpLunVwdP1Cv382rlFdToTFvXgP28vdy1dG64UpMjtGh26KBjYIbD\nFTMPA8BYQZgBvobm1i7tO1qrvYU1OnC8bsC1kKSzK1Ivmxeu5cmRumLmVHm4j+xTSM6aeRgAxiLC\nDDBMVfUt2lNYoz2FNTpa0qA++8D9pgRMUuq8CC1PjtTcuGCZzaO3IvVozzwMAGMZYQa4iN4+u46V\nntLewhrtPVIja93AywhIUsQUX6UmRWh5coRmRgfJzYkz8Y70zMMAYBSEGWAA7Z09+ux4nfYU1mjf\n0Vo1tw48/sVkkmbHBGlpYrhSEsMVHebHWkgA4GSEGeBLdafatO/Y2fEvh07a1NPbN2A/Tw+zFswK\nUUpiuBbPDVOQ3yQnVwoAOB9hBhNWd0+fjpY2aP/ROu07VjvoBHbS2ZWolyaGa2liuK6YGcIyAgAw\nhhBmMKE0NLVr/7E67Ttaq4Mn6gd9+kiSpkX4K+XLADPDEujU8S8AgEtHmMG41tvbp2Nljdp/rFb7\njtaqpGrgFaglycPdTUnxU7VoTqhSEiMUFuzjxEoBAJfLJWHmmWee0aFDh2QymbRx40YlJSU52vLy\n8vTcc8/JbDYrLi5OmzdvdkWJMDDb6XYdPFGn/cfq9NmJerW2dw/aNzTIW4vnhGnRnDAlx08dkQns\nAADO5fTv3Pn5+SorK1NmZqaKioq0adMmZWZmOtqfeOIJvfrqqwoNDdWDDz6ojz/+WFdddZWzy4SB\ndHT2qKC4QZ+dqNNnx+tVUTv42Bd3s0mJ06ecDTAJYbKETubpIwAwOKeHmdzcXGVkZEiS4uPj1dzc\nrNbWVvn6+kqStm7d6vhzcHCwTp8+7ewSMcb19dlVXNmkz07U6eCJeh0pOTXok0eSNDVgkhbNCdPi\nOWFKnjH1oitQAwCMxelhxmazad68eY7toKAg2Ww2R4A593tdXZ1ycnL00EMPObtEjEH1jWdvHR08\nUa+DJ+sHnfdFOjv2JTFuiubPCtGiOWGKDWfuFwAYz1w+QMBu7z8XfENDg+655x49+eSTCggIcEFV\ncLXm1i4VFNl0+AubDp2sH3LWXensk0fzZ4VowaxQzZ0erEmeLv+nDQBwEqd/xw8NDZXNZnNs19XV\nKSQkxLHd0tKiH/zgB3r44YeVmprq7PLgIm0d3SosbtDhL2w6fNKmkuomDZBzHQIne50NL7NDdMXM\nEE0J8HZesQCAMcXpYSYtLU0vvfSSbr31VhUWFiosLEw+Pv98BHbLli268847lZaW5uzS4ESd3b06\nVnJKh76o1+EvbDpZcVp9g63YqH/eOlowO0TzZ4VqWoQ/874AACS5IMwsWLBAiYmJWrt2rcxmsx5/\n/HFt27ZNfn5+WrFihf72t7+pvLxcb775pkwmk775zW/qlltucXaZGGGd3b06Ud6ogqIGHf6iXsdK\nG4cctOvmZtLM6EAlz5iq5BlTlTCNW0cAgIGZ7AMNWjEQq9Wq9PR07dq1SxaLxdXl4EttHd06WnpK\nhcUNKixu0Iny00OGF5NJiosIUPLMs+ElcfoUnjoCAFwS/quLEdHU0qkjJQ0qKG7QkeIGFVc2aYi7\nRpKk6LDJSp4RouQZUzUvfqr8fT2dUywAYFwhzOCy2E63q+DLqy6FxQ1DTlR3jiV0shKnT1FS/FQl\nzZiqYH9WmwYAfH2EGVxUb2+fymrO6GjpKR0rO6WjJadUe6ptyNeYTFJcZIASp09R4vQpmhsXrCA/\nwgsAYOQRZtDPmbYuHS9rPBteSk/pRHmjOrp6h3yNu9mkmdFBmhsXrHnxZwfsTvZmzAsAYPQRZia4\nvj67KurO6Fhpo459eeXlYhPUSZKnh1lzpgUpMW6KEuOnaFZMEE8bAQBcgp8+E8yZti6dLD+t42Wn\ndKysUcfLTqm1o+eir5sSMEkJ04I1Z1qwEmKDND0qUB7ubk6oGACAoRFmxrGOzh4VVTbpZEWjTpaf\n1smK06puaL3o68xuJsVbApQQG6yEacFKiA1WSBAz7AIAxibCzDjR09unsupmnag4rZPljTpZcVrl\nNc0XfTxaOrs0QMK0IEd4mREdKC8P8+gXDQDACCDMGFBvn11V9S36wnr2asuJ8kaVVDapq2fwSenO\ncTebNC0yQDOjA7+8ZRSs8Ck+rCoNADAswswY193Tp/KaZhVXNqmosknFlU0qqWq66NNF0tnHoy2h\nfpoZHahZMUGaGR2ouEh/ebhz1QUAMH4QZsaQjs4elVQ1q7jytIq+DC/lNc3q6b20FSdCg7w1MzpI\ns2ICNTM6SPGWAJYEAACMe4QZF2k806Gy6guvuFTWt+hSV8oK8vNSvCVQM6PP/QpSoJ/X6BYNAMAY\nRJgZZZ3dvaqoOaPS6iaVVp/9vaz6jE63dF7ye4QF+yjeEqDpUQGKjwrU9KgAlgIAAOBLhJkR0tdn\nV11jm0qrm//5q6pZ1baWS3qiSJLcTJIlzO/L0HI2uMRFBTCTLgAAQyDMDJPdbldDU4fKa8+o4stf\npdXNKq9pVnvnxQflnuPlada0cH9Ni/R3hJfYCH9m0QUAYJj4yTmIvj676k+3q6L2jMpr/hlcymvP\nqL3z4jPmnmMySZFTfRUb4a9pEQGaFuGnaREBCgv2kZsbj0MDAPB1Tfgw09vbp9rGNllrWxxXW8pr\nz8hae+aSHn8+n7+vp6ZF+Dt+xUb4Kybcj6stAACMognxU9Zut+t0S6cq61pUWd+qqvoWVX75q6ah\n9ZIffT7H19tDMWF+iv7yV0y4n+Ii/BXo58XkcwAAONm4CjPtnT2qqm9RVX2rrPUtqqpvcfzedgmL\nKX6Vv6/n2bByLrSE+Sk63E9BhBYAAMaMcRNmHvmvj9Xa63tZr50SMElRIZMVE+53wRWXgMnM2wIA\nwFg3bsLM6TOd8vAZPMx4e7krKnSyLCGTFRkyWVEhvor68s/eXuPmrwEAgAlnXP0UN7uZFD7FV5bQ\nc4Hln6GF8SwAAIxP4ybMPH1vmq6YO0Nms5urSwEAAE40bn7yhwb5EGQAAJiA+OkPAAAMjTADAAAM\njTADAAAMjTADAAAMjTADAAAMjTADAAAMjTADAAAMjTADAAAMjTADAAAMjTADAAAMjTADAAAMjTAD\nAAAMjTADAAAMjTADAAAMzd0VB33mmWd06NAhmUwmbdy4UUlJSY62nJwcPffcczKbzbrqqqt07733\nuqJEAABgEE6/MpOfn6+ysjJlZmbqqaee0ubNmy9o37x5s1566SW98cYbys7OVlFRkbNLBAAABuL0\nKzO5ubnKyMiQJMXHx6u5uVmtra3y9fVVRUWFAgMDFRYWJklauXKl8vLyFB8f7+wyMQqyy/O17UiW\nrM3VsvhHaM3c1UqLWeLqsgAABuf0MGOz2TRv3jzHdlBQkGw2m3x9fWWz2RQcHOxoCw4OVkVFhbNL\nxCjILs/XC7mv6NNtn6qjqlmTIv1VvqZSksZloFm/fr1KSkoUFxen3/zmN0477kMPPaTi4mJNnz5d\nzz//vNOOC+MzmUyOP9vt9kH3DdXfGVx1XCO56667dPz4cc2ePVuvvPKKq8txCpcPAB7qHyP/UMeP\nbUey9Om2T+XXY1ZIaJD8esz6dNunevtIlqtLG3Hr169Xe3u7wsPD1d7ervXr1zvluA899JDa2toU\nHh6utrY2PfTQQ045Lozv/IBwbnugfUP1dwZXHddI7rrrLrW1tSk6OlptbW266667XF2SUzg9zISG\nhspmszm26+rqFBIS4mirr693tNXW1io0NNTZJWIUWJur1VHVfMG+jqpmWZurXVTR6CkpKRlye7QU\nFxcPuQ1g/Dt+/PiQ2+OV08NMWlqasrLO/m+8sLBQYWFh8vHxkSRFRUWptbVVVVVV6unp0UcffaQV\nK1Y4u0SMAot/hCZF+l+wb1Kkvyz+ES6qaPTExcUNuT1apk+fPuQ2gPFv9uzZQ26PV04PMwsWLFBi\nYqLWrl2rp59+Wo8//ri2bdumnTt3SpKeeOIJrV+/Xt/97nd10003KTY21tklYhSsmbtaK9as0Bn3\nXtXXNeqMe69WrFmhm+eudnVpI+43v/mNvL29VVNTI29vb6eNmXn++efl4+Ojmpoa+fj4MGYGl+yr\nt/TtdvuA+4bq7wyuOq6RvPLKK/Lx8VFFRYV8fHwmzJgZk93g/xqsVqvS09O1a9cuWSwWV5eDIWSX\n5+vt855mupmnmQAAI8Alk+ZhYkqLWUJ4AQCMOJc/zQQAAPB1EGYAAIChEWYAAIChEWYAAIChEWYA\nAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAICh\nEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYA\nAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAICh\nEWYAAIChEWYAAIChEWYAAIChEWYAAIChEWYAAIChuTv7gD09PXrsscdUVVUls9msZ555RhaL5YI+\nf//73/U///M/MpvNSklJ0Y9//GNnlwkAAAzC6Vdm3n33XQUEBOj111/XD3/4Qz377LMXtHd0dOjZ\nZ5/Vn//8Z2VmZio3N1dFRUXOLhMAABiE08NMbm6uMjIyJEnLly/XgQMHLmifNGmS3nnnHXl7e0uS\nAgMDdfr0aWeXCQAADMLpYcZmsyk4OFiSZDKZ5Obmpp6engv6+Pj4SJKOHz+uqqoqzZ8/39llAgAA\ngxjVMTN/+ctf9Ne//lUmk0mSZLfbdfjw4Qv69PX1Dfja0tJSPfLII3r22WdlNpsHPUZvb68kqaam\nZoSqBgAArhYeHi5390uLKaMaZm655RbdcsstF+zbsGGDbDabZs+e7bgi89Via2pq9MADD+jXv/61\nZs+ePeQx6uvrJUm33XbbCFYOAABcadeuXf0eEBqM059mSktL044dO5SWlqYPPvhAKSkp/fps2rRJ\nTzzxhBISEi76fvPmzdNrr72mkJCQIa/gAAAA4wgPD7/kvia73W4fxVr66evr06ZNm1RWViYvLy9t\n2bJFYWFh+v3vf6+UlBQFBARozZo1SkpKkt1ul8lk0p133qmrr77amWUCAACDcHqYAQAAGEnMAAwA\nAAyNMAMAAAyNMAMAAAzN6U8zfV0dHR167LHH1NDQoK6uLt1zzz1KSEjQT37yE9ntdoWEhOhXv/qV\nPDw8XF0qRshAn3lWVpYKCgoUFBQkSbr77ru1cuVKF1eK0dDZ2ambbrpJ9913n5YtW8a5PgGc/5nv\n2bOHc32c27t3rx588EHNnDlTdrtds2fP1ve///1hneuGCzMffPCBkpKSdPfdd6uqqkp33nmnFi5c\nqO9+97tavXq1nnvuOb311ltau3atq0vFCBnsM3/kkUf4pjYBvPzyywoMDJQkvfDCC7r99tt17bXX\ncq6PY+d/5pI41yeApUuX6oUXXnBsb9iwYVjnuuFuM91www26++67JUlVVVWKiIhQfn6+Vq1aJUm6\n+uqrlZOT48oSMcIG+sylszNKY3wrLi5WcXGxVq5cKbvdrvz8fMc0DZzr49NXP3OJc30i+OpnvHfv\n3mGd64YLM+esXbtWjz76qDZs2KD29nbH5acpU6Y4ZgXG+HLuM9+4caMk6bXXXtMdd9yhhx9+mMVI\nx6lf/vKXeuyxxxzbnOvj3/mf+bmlcDjXx7+ioiLde++9uu2225STk6OOjo5hneuGu810TmZmpo4d\nO6ZHHnnkgkRHgh+/zv/MN27cqMDAQCUkJOj3v/+9XnzxRf3sZz9zdYkYQW+//bYWLFigqKioAds5\n18efr37mdrtd3/rWtzjXx7nY2Fjdf//9uv7661VRUaHvfe97FyxAfSnnuuHCTGFhoaZMmaLw8HAl\nJCSor69Pvr6+6urqkqenp2praxUaGurqMjGCvvqZ9/b2atasWY7V19PT0/Xkk0+6tkiMuN27d8tq\nterDDz9UbW2tPDw85OPjw7k+jp3/mdfU1MjLy0s///nPHUvbcK6PT2FhYbr++uslSdHR0Zo6daoK\nCgqGda4b7jZTfn6+XnnlFUmSzWZTW1ubUlNTtWPHDklSVlaWrrzySleWiBE20Gf+xBNPqKKiQpK0\nZ88ezZo1y5UlYhQ899xz+stf/qL/+7//07/+67/qvvvu41wf587/zG+55Rbde++9euONNzjXx7l3\n3nnH8T2+vr5eDQ0N+va3vz2sc91wyxl0dnZq48aNqqmpUWdnpx544AElJibq0UcfVVdXlyIjI/XM\nM8+w6OQ48tXP/P7775ePj49+9atfydvbW76+vnr66acdV2ow/rz00kuyWCxasWIF5/oEce4zj4yM\n5Fwf51pbW/Xwww/rzJkz6unp0f3336+EhAT99Kc/veRz3XBhBgAA4HyGu80EAABwPsIMAAAwNMIM\nAAAwNMIMAAAwNMIMAAAwNMIMAAAwNMIMgFFVWVk5rBWPd+/erebmZknSww8/rLq6OknS3/72t2Ed\n9/nnn9dLL700rNcAMCbCDIBRd27BwEvxv//7v47FBJ999lmFhoaqtrZWmZmZo1UeAIMjzABwif/6\nr//SunXrtG7dOq1fv169vb164403tG/fPv3kJz9RUVGRVq1apYqKCj3yyCM6efKkHnvsMe3du1fr\n1q1zvM+GDRv017/+VdLZ6fC/9a1v6Yc//KHKysocffLy8nT77bfr9ttv11133aXKykqnf70ARg9h\nBoDT9fb2ytvbW6+//rpef/11NTU16dNPP9V3vvMdTZ06Vf/5n/+p+Ph4xxWdBx54QLNmzdKWLVsk\nDXylp7S0VO+++67eeust/fa3v3WEmY6ODj355JP67W9/q1dffVW33XabfvnLXzrviwUw6gy3ajYA\n4zObzXJzc9N3vvMdubu7q6SkRI2Njf36DWe1lRMnTigxMVHu7me/rS1evNixv76+Xvfff7/sdrvs\ndvuwbnsBGPsIMwCc7sCBA9q6dau2bt0qLy8v/ehHP7rk1341iHR1dUk6G3zc3P55sbmvr0+S5Onp\nqcjISP35z38egcoBjEXcZgIw6r56haWhoUFRUVHy8vJSZWWlDh065Aglbm5u6u7uvqD/+fsmT56s\n2tpaSVJ7e7sOHz4sSYqPj9eRI0fU09Oj7u5u7d27V5IUFxenxsZGnTx5UpKUn5+vN998c/S+WABO\nx5UZAKOusbFR3/ve9xy3eBITE9XS0qLbbrtNM2bM0AMPPKCXX35ZKSkpWrFihe655x5t2bLFcRVm\nxowZstlsuvvuu/XHP/5Rs2bN0re//W3FxMRo4cKFjj7p6em69dZbFRkZqblz50qSvLy89Otf/1qb\nNm2Sl5eXJOkXv/iFa/4iAIwKk304N6UBAADGGG4zAQAAQyPMAAAAQyPMAAAAQyPMAAAAQyPMAAAA\nQyPMAAAAQyPMAAAAQyPMAAAAQ/v/s6Jybuepo4YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8812e3bba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plot_log_reg(x='Latitude', y='Allele', data=raw_df, clf=clf, xmin=30, xmax=50, alpha=0.02)\n",
    "df.sort_values('Latitude').plot('Latitude', 'p, Mpi100', ls='', marker='o', ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.19685284] [[ 0.]]\n"
     ]
    }
   ],
   "source": [
    "clf0 = log_reg_null_model(raw_df['Allele'])\n",
    "print(clf0.intercept_, clf0.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.0572187317602719e-20"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_reg_lik_ratio_test(raw_df['Latitude'], raw_df['Allele'], clf0, clf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}